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Theories  of  motion  and  stability  of  fluid-saturated  soils,  including  the  commonly 
used  engineering  approach  to  liquefaction  analysis  as  well  as  theories  based  on  mechan¬ 
ics  of  mixtures,  are  critically  examined.  Description  erf  motion,  development  of  equations 
of  balance,  constitutive  relationships  as  well  as  development  of  solution  procedures  are 
reviewed.  Limitations  of  various  theories,  their  similarities  as  well  as  inconsistencies  are 
identified.  Laboratory  investigations  into  dynamic  behavior  of  saturated  soils  are 
reviewed. 
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A  theory  of  dynamics  of  saturated  soils  using  a  convected  coordinate  system  to 
describe  the  motion  of  9oil  particles,  and  describing  the  motion  of  the  fluid  as  relative 
to  the  solid,  is  described.  '  This  theory/-  may  be  regarded  asan  extension  of  Gibson's 
theory  of  non-linear  soil  consolidation  to  three-dimensions  and  to  include  inertia  effects. 

Solution  procedures ,  developed  for  certain  specializations  of  the  equations  of  motion 
of  saturated  soils,  are  described.  These-  include  analytical,  semi-analytical  and  numerical 
solution  schemes.  The  finite  element  is  selected  as  the  numerical  procedure  for  approxi¬ 
mate  solution.  Spatial  discretization,  time  domain  solution  procedures  and  alternative  for¬ 
mulations  of  the  field  equations  through  a  variational  formulation  are  discussed. 

Shaking  table  tests  for  validation  of  various  theoretical  concepts,  performed  on 
saturated  Ottawa  sand'  are  described.  These  included  tests  on  anisotropically  as  well  as 
isotropically  consolidated  samples  and  tests  to  study  the  effect  of  overburden  on  a  soil 
system  subjected  to  shaking.  Harmonic  as  well  as  frequency  banded  random  amplitude 
excitations  were  used.  Results  are  compared  with  previous  laboratory  investigations. 


TABLE  OF  CONTENTS 


FOREWORD .  i 

ABSTRACT  .  iii 

Section  L  INTRODUCTION  .  1 

THE  PROBLEM  .  1 

RESEARCH  OBJECTIVES  . 2 

RESEARCH  PLAN.  . 4 

ANALYTICAL  STUDIES . 4 

LABORATORY  INVESTIGATIONS . 5 

COMMENTS . 5 

STRUCTURE  OF  THE  REPORT  .  8 

Section  H:  REVIEW  OF  PREVIOUS  WORK .  11 

INTRODUCTION .  11 

THE  ENGINEERING  APPROACH .  .  11 

CONTINUUM  MECHANICS  APPROACHES .  13 

Biots  Theory .  13 

Theories  of  Mixtures. .  15 

Introduction  .  15 

Density  of  Each  Constituent  and  of  the  Mixture.  .  16 

Description  of  Motion. .  16 

Measures  of  Deformation. .  18 

Balance  Laws .  20 

Constitutive  Relations .  34 

Comments. .  51 

Section  HL  A  DYNAMICAL  THEORY  OF  SATURATED  SOILS .  59 

INTRODUCTION .  59 

KINEMATICS .  60 

BALANCE  LAWS .  61 

Mass  Balance  of  the  Solid .  61 

Mass  Balance  of  the  Fluid .  62 

Balance  of  Momentum  of  the  Fluid  Phase .  64 

Balance  of  Momentum  of  the  Fluid-Saturated  Solid .  67 

Balance  of  Angular  Momentum  of  the  Fluid  Saturated  Solid .  69 

SOME  SPECIALIZATIONS .  70 

Specialization  to  One-Dimensional  Problem .  70 

Small  deformation  Theory .  73 


CONSTITUTIVE  RELATIONS 


74 


Section  IV:  SOLUTION  PROCEDURES  .  77 

INTRODUCTION .  77 

EXACT  SOLUTIONS  . 78 

Introduction .  78 

Craig's  Solution .  79 

Integration  of  Craig's  Solution .  86 

Evaluation  of  Garg's  Approximations  .  88 

Wave  Propagation  in  a  Fluid-Saturated  Soil  Layer .  93 

SEMI-DISCRETE  SOLUTION . 94 

FINITE  ELEMENT  SOLUTIONS .  97 

Variational  Formulation  .  97 

Two-  and  Three-Field  Formulations . 103 

Spatial  Discretization  . 106 

Singularity  Elements  . 108 

Time  Domain  Integration . 108 

Nonlinear  Problems  . 133 


Section  V:  LABORATORY  INVESTIGATIONS . 145 

INTRODUCTION . 145 

EXPERIMENTAL  FACILITIES . 148 

Introduction . 148 

.Shaking  table . 148 

Test  Box . 149 

Rubber  Membrane  . 152 

Instrumentation . 152 

Test  Material . 152 

TESTING  PROCEDURES . 155 

Sample  saturation  . 155 

Sample  Confinement . 155 

Input  Motions . 156 

Harmonic  motion . 156 

Random  Motion . 158 

RESULTS  . 158 

DISCUSSION  . 165 

Non-uniformity  of  sand  samples . 167 

Nonuniformities  due  to  testing  . 167 

Membrane  penetration . 168 

Nonuniformity  of  the  confining  pressure  at  the  sample  boundaries  . 169 


Section  VL*  DISCUSSION . 171 

OBJECTIVE  OF  THE  RESEARCH  PROGRAM . 171 

ACCOMPLISHMENTS . 171 

Review  of  Theories  . 171 

A  Dynamical  Theory  of  Saturated  Soils. . 174 

Development  of  Solution  Procedures. . 174 


vi 


Laboratory  Investigations. . 176 

FUTURE  WORK . 177 

Section  VIL  LIST  OF  REFERENCES . 179 

Appendix  A:  CONVECTED  COORDINATES  . 191 

NOTATION . 191 

COORDINATE  TRANSFORMATION  . 191 

CONTRA  VARIANT  AND  COVARIANT  VECTORS . 193 

CURVILINEAR  COORDINATES . .194 

BASE  VECTORS  AND  METRIC  TENSORS . 194 

PHYSICAL  MEANING  OF  COVARIANT  AND  CONTRA  VARIANT 

COMPONENTS  OF  A  VECTOR . 198 

PARTIAL  DERIVATIVES  IN  CARTESIAN  COORDINATES  . 201 

COVARIANT  DIFFERENTIATION  . 202 

GEOMETRIC  INTERPRETATION  OF  COVARIANT  DERIVATIVES . 203 

GRADIENT.  DIVERGENCE  AND  CURL  IN  CURVILINEAR 

COORDINATES . 204 

KINEMATICS  AND  KINETICS  . 206 

Introductory . 206 

Geometric  Relations  . 206 

Description  of  Deformation  . 209 

Velocity  and  Acceleration  in  Convected  Coordinates . 214 

Changes  in  Volume  and  Area  . 215 

Kinetics  . 216 

Appendix  B:  PUBLICATIONS  AND  PRESENTATIONS . 221 

RESEARCH  REPORTS  . 221 

CONFERENCE  PROCEEDINGS . 223 

REFEREED  JOURNAL  ARTICLES . 223 

DISSERTATIONS  AND  THESES  . 224 

OTHER  PRESENTATIONS  . 224 


LIST  OF  TABLES 


Mechanical  Properties  Und  in  the  Example 


LIST  OF  FIGURES 


Variable*  Involved  in  Dynamics  of  Saturated  Soil*  . 

Research  Accomplished  Under  AFOSR  Grant  83-0055 . 

Geometry  of  an  Infinitesimal  Parallelepiped . 

Fluid  Equilibrium  in  the  Deformed  State  . 

Equilibrium  of  the  Fluid-Saturated  Reference  Volume  in  the  Current 
Configuration  . 

Velocity  Excitations  Applied  at  the  Boundary  . 

Family  of  Variational  Principles  by  Direct  Formulation  . 

Family  of  Complementary  Variational  Principles . 

Representative  Sail  Column  . 

Types  of  Excitation  Applied  on  the  Soil  Column  . 


Example 

1: 

Solid 

Velocity 

History 

at 

10 

cm. 

(Low 

Drag) 

Example 

1: 

Fluid 

Velocity 

History 

at 

10 

cm. 

(Low 

Drag) 

Example 

1: 

Solid 

Velocity 

History 

at 

30 

cm. 

(Low 

Drag) 

Example 

1: 

Fluid 

Velocity 

History 

at 

30 

cm. 

(Low 

Drag) 

Example 

1: 

Solid 

Velocity 

History 

at 

10 

cm. 

(High 

Drag) 

Example 

1: 

Fluid 

Velocity 

History 

at 

10 

cm. 

(High 

Drag) 

Example 

1: 

Solid 

Velocity 

History 

at 

30 

cm. 

(High 

Drag) 

Example 

1: 

Fluid 

Velocity 

History 

at 

30 

cm. 

(High 

Drag) 

Example 

2: 

Solid 

Velocity 

History 

at 

10 

cm. 

(Low 

Drag) 

Example 

2: 

Fluid 

Velocity 

History 

at 

10 

cm. 

(Low 

Drag) 

Example 

2: 

Solid 

Velocity 

History 

at 

30 

cm. 

(Low 

Drag) 

Example 

2: 

Fluid 

Velocity 

History 

at 

30 

cm. 

(Low 

Drag) 

Example 

2: 

Solid 

Velocity 

History 

at 

10 

cm. 

(High 

Drag) 

Example 

2: 

Fluid 

Velocity 

History 

at 

10 

cm. 

(High 

Drag) 

Example 

2: 

Solid 

Velocity 

History 

at 

30 

cm. 

(High 

Drag) 

26  Example  2:  Fluid  Velocity  History  at  30  cm.  (High  Drag)  . 132 

27.  Discretization  of  the  Soil  Column  . 136 

28.  Sress  History  for  Case  3.  . 137 

29.  Solid  Column  Under  Stress  Pula  of  Short  Duration  . 138 

30.  Strea  Profile  at  t  -  4. . 139 

31.  Strea  Profile  at  t  -  8. . 140 

32.  Strea  Profile  at  t  -  12. . 141 

33.  Strea  Profile  at  t  -  16. . 142 

34.  Strea  Profile  at  t  -  20 . 143 

35.  Large  Scale  Liquefaction  Testing  System . 150 

36.  Sample  Test  Chamber . 151 

37.  Relative  Positions  of  Pore  Preanin  Transducers . 153 

38.  Ottawa  Sand  Gradation . 154 

39.  Comparison  of  Table  and  Sample  Acceleration  Time  Histones  for 

Harmonic  Excitation . 157 

40.  Fourier  Spectrum  at  Random  Amplitude  Acceleration . 159 

41.  Liquefaction  Test  Results  for  Harmonic  Excitation . 161 

42.  Liquefaction  Test  Results  for  Anisotropically  Consolidated  Samples  . 162 

43  Liquefaction  Test  Results  Obtained  with  a  Reaction  Mass  on  the 

Sample  . 163 

44.  Liquefaction  Test  Results  for  Random  Amplitude  Excitation . 164 

45.  Coordinate  Surfaces  . 195 

46  Position  Vector  and  Its  Differential . 196 

47.  Covariant  and  Contra  variant  Components  of  a  Vector  . 200 

48.  Deformation  of  a  Body . 207 

49.  Converted  Coordinate  System  . 210 

50.  Infinitesimal  Curvilinear  Tetrahedron . 217 


-  xu  • 


Section  I 


INTRODUCTION 


1.1  THE  PROBLEM 

Dynamic  loading  of  fluid-saturated  geological  deposits  result!  in  changes  in  the 
fluid  pressure  as  well  as  the  stress  field  in  the  solid  matrix.  This  phenomenon  can.  in 
some  cases,  lead  to  instability  of  the  soil  matrix  resulting  in  ami  liquefaction.  Addition¬ 
ally,  the  energy  dissipation  associated  with  the  relative  oscillatory  motion  between  the 
fluid  and  the  solid  could  introduce  attenuation  of  the  propagating  wave.  This  effect 
would  be  in  addition  to  the  energy  dissipation  due  to  any  inelastic  deformation  of  the 
soil  matrix.  Inasmuch  as  the  tranamissibility  of  motion  through  the  soil  layer  depends 
upon  the  soil  characteristics,  the  soil  deposit  may  act  as  a  selective  filter/amplifier. 
For  prediction  of  behavior  of  m curated  soils  under  dynamic  loads  due  to  blast,  earth¬ 
quake  or  other  dynamic  event,  it  is  important  to  develop  adequate  methodology. 

Reliable  analysis  of  saturated  soil  deposits  subjected  to  dynamic  loads  involves  the 
following  three  steps: 

1.  Correct  formulation  of  the  equations  of  dynamic  equilibrium. 

2.  Correct  representation  of  material  behavior. 

3.  Exact  or  approximate  solution  of  the  problem. 

The  theoretical  model  must  be  validated  by  laboratory  and/or  field  testa.  As  exact 
solutions  to  the  boundary  value  problems  represented  by  the  laboratory  and  field  tests 
may  not  be  available,  it  may  be  necessary  to  use  numerical  solution  procedures.  How¬ 
ever,  before  any  computational  technique  can  be  used  with  confidence,  it  must  be  veri- 
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tied  against  exact  solutions  to  the  theoretical  model.  It  may  be  necessary  to  construct 
exact  solutions  to  tome  simple  problems  for  code  verification. 

In  view  of  the  apparent  diversity  of  opinions,  postulates,  and  assumptions  made 
in  setting  up  various  approaches  to  the  problem,  there  was  evident  need  to  develop 
more  realistic  models  of  mechanical  behavior  of  fluid-saturated  soils.  The  research  pro¬ 
gram  supported  by  the  Air  Force  Office  of  Scientific  Research  was  designed  to  address 
all  the  three  components  of  the  problem. 

lJ  RESEARCH  OBJECTIVES 

Figure  1  illustrates  the  variables  involved  in  the  problem  of  dynamics  of  fluid- 
saturated  soils.  In  the  absence  of  water,  the  dry  soil  would  be  modeled  as  a  contin¬ 
uum.  Its  motion  would  be  defined  completely  by  the  components  of  displacement, 
velocity  and  acceleration  designated  u,,  u,,  u,,  in  the  figure.  These  would  define  strains, 
and  strain-rates,  which  would  be  related  through  constitutive  Laws  to  the  stresses,  and 
the  stress  rates.  Equilibrium  relations  would  relate  the  stresses  with  the  applied  body 
and  surface  forces.  Similarly  the  mechanics  of  the  fluid,  in  isolation  from  the  solid  are 
well-known.  However,  when  the  two  occur  together,  the  problem  immediately  becomes 
much  more  complex.  The  equilibrium  equations  for  each  of  the  two  constituents  may 
involve  interaction  terms;  the  streaes  in  each  buy  depend  upon  the  kinematics  of  both; 
and  there  may  be  other  couplings  in  the  behavior.  Various  theories  have  been  proposed 
from  time  to  ume  to  explain  the  mechanical  behavior  of  saturated  soils  and  methodolo¬ 
gies  have  been  suggested  for  analysis  of  liquefaction. 

The  objective  of  the  research  program,  waa  to  critically  evaluate  the  current 
methods  of  analysis  in  the  light  of  recent  developments  in  theories  of  interacting  con¬ 
tinue  and  to  develop  thermodynamically  and  physically  consistent  theoretical  or  mathe- 


A.  inertial  coupling 

B.  CONSTITUTIVE  COUPLING 

C.  BALANCE  COUPLING 

(MASS. MOMENTUM. ENERGY) 


Figure  1:  Variables  Involved  in  Dynamics  of  Saturated  Soils 
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models  of  fluid-saturated  soil  systems.  Theories  of  interacting  continue,  in  gen¬ 
eral.  regard  a  fluid-saturated  soil  as  a  multi -component  mixture  (superposed  continual. 
This  approach  had  been  successfully  implemented  in  the  analysis  of  static  and  quasi¬ 
static  response  of  saturated  soils.  Several  investigators  have  tried  to  extend  this 
approach  to  the  dynamic  case.  However,  there  have  been  deficiencies  in  realistic  nmula- 
tion  of  behavior  of  saturated  toils.  The  purpose  of  the  research  program  waa  to  review 
the  theoretical  basis  of  equations  governing  behavior  of  soil-water  mixtures  under 
dynamic  conditions,  including  interaction  between  soil  and  water.  The  theoretical 
development  was  to  be  implemented  in  effective  finite  element  computer  programs 
incorporating  recent  developments  in  coding  to  ensure  optimal  combination  of  solution 
accuracy  and  economy.  The  analytical  research  waa  to  be  supplemented  by  a  program 
of  laboratory  investigations. 

lJ  RESEARCH  PLAN. 

A  two-year  program  of  research  was  initially  approved  by  the  AFOSR  starting 
February  1,  1983.  This  was  later  extended  to  a  four  year  effort  ending  January  31, 
1987  and  finally  by  another  one  year  and  one-month  to  February  29,  1988.  The 
research  plan  included  both  theoretical  and  laboratory  investigations. 

\A  ANALYTICAL  STUDIES 

The  first  step  in  the  research  program  was  to  carefully  examine  the  theoretical 
underpinnings  of  various  existing  theories  of  motion  and  stability  of  fluid-saturated 
soils.  This  investigation  covered  a  range  of  theories  and  procedures,  including  the  com¬ 
monly  used  engineering  approach  to  liquefaction  analysis  and  extending  to  mathematical 
theories  based  on  mechanics  of  mixtures.  A  theory  properly  describing  the  physical  phe¬ 
nomena  was  to  be  selected/developed  and  implemented  in  a  solution  procedure  to  pre¬ 
dict  the  liquefaction  behavior  of  the  saturated  soil  under  dynamic  loading. 
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1.3  LABORATORY  INVESTIGATIONS 

To  provide  a  database  for  validation  of  various  theoretical  concepts,  shaking  table 
tests  were  to  be  performed  on  saturated  Ottawa  sand.  A  program  of  tests  on  anisotrop- 
icaliy  as  well  as  isotropically  consolidated  samples  was  to  be  carried  out. 

1  j6  COMMENTS 

The  research  program  essentially  followed  an  evolutionary  approach.  As  the  first 
step  in  the  program,  literature  on  the  subject  was  carefully  reviewed.  Two  approaches 
were  selected  for  detailed  investigation.  These  were  the  popular  "engineering  approach" 
introduced  by  Seed  and  his  co-workers  [77-80,148]  for  liquefaction  studies,  and  the 
theories  of  machine*  of  mixtures  including  Biot's  theory  [17-19]  which  has  been  the 
basis  of  analytic  and  numerical  solutions  for  the  problem.  A  dynamical  theory  for 
dynamics  of  saturated  soil  was  developed  as  an  extension  of  Gibson's  theory  of  non¬ 
linear  consolidation.  For  interpretation  of  test  data  in  terms  of  predictions  from  vari¬ 
ous  theories,  it  was  necessary  that  solutions  to  boundary  value  problems  defined  by  the 
competing  theories  be  available.  It  was  found  that  solutions  to  only  the  simplest  con¬ 
figurations  were  available  for  Biot's  theory.  Approximate  solution  schemes  had  been  pro¬ 
posed  by  several  investigators  but  these  had  not  been  adequately  verified.  Finite  ele¬ 
ment  computer  codes  were  developed  for  analysis  of  dynamic  response  of  saturated  soils 
for  linear  as  well  as  nonlinear  material  properties  for  the  engineering  approach  as  well 
as  Biot's  theory  and  its  appropriate  extension.  Semi-analytic  solution  procedures  were 
developed  to  serve  as  bench-marks  for  validation  of  the  time-domain  integration 
schemes.  Analytical  solutions  were  developed  for  certain  problems  for  the  purpose  of 
code  verification. 

Validation  of  numerical  procedures  presented  a  difficulty  because  of  the  paucity 
of  exact  solutions  even  for  relatively  simple  problems.  Garg's  [50]  fundamental  solution 
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for  Biot's  equations  specialized  one-dimensional  wave  propagation  in  fluid-saturated 
media  was  integrated  to  develop  solutions  for  several  cases  of  dynamic  loading  on  the 
surface  of  a  saturated  soil  cnim""  Inconsistencies  in  the  solution  for  'strong'  coupling, 
for  which  the  drag  between  the  soil  and  the  fluid  is  high,  were  removed  but  it  was 
noticed  that,  in  order  to  obtain  an  exact  solution  to  the  problem,  Garg  had  made 
assumptions  which  reduced  the  value  of  his  theoretical  solutions  u  bench-marks  to  ver¬ 
ify  numerical  procedures.  To  meet  the  needs  of  the  research  program  effectively,  correct 
solutions  to  this  problem  were  developed.  It  was  found  that  for  the  materials  com¬ 
monly  encountered  and  for  a  short  period  of  time  after  application  of  a  sudden  (e.g. 
blast)  load,  Garg's  approximation  is  quite  good.  Analytical  solutions  for  some  simple 
cases  were  developed  for  propagation  of  standing  waves  by  separating  the  singularity 
from  the  smooth  diffusive  wave  motion  and  using  a  combination  of  the  method  of 
characteristics  and  finite  element/finite  difference  procedures.  These  solutions  were 
extended  to  some  two-dimensional  cases.  Computer  codes  were  tested  against  these  exact 
solutions. 

To  develop  alternative  finite  element  approaches,  a  variational  formulation  of 
Biot's  theory,  along  with  various  extensions  and  specializations,  was  developed.  It  served 
as  the  basis  for  the  two-field  and  the  three-field  finite  element  solution  procedures. 
Elements  suitable  for  wave  propagation  analysis  were  selected/developed.  Singularity 
elements  had  to  be  used  near  loaded  boundaries.  For  nonlinear  problems,  an  incremen¬ 
tal  approach  was  necessary.  The  equations  governing  this  case  along  with  a  variational 
formulation  of  the  problem  were  developed.  Only  material  nonlinearity  was  considered. 
The  theory  was  implemented  in  a  finite  element  computer  program.  A  modular  struc¬ 
ture  was  used  so  that  a  variety  of  models  could  be  selected.  The  well-established  'cap 
model'  was  implemented  to  fit  the  behavior  of  the  sand  used  in  the  laboratory  experi¬ 
ments.  For  constitutive  equations  in  incremental  form,  the  corresponding  balance  equa- 


tinn«  too  were  written  in  terms  of  incremental  quantities.  This  introduced  incremental 
changes  in  porosity  and  mass  density  u  additional  variables.  Nonlinear  codes  were 
checked  against  available  solutions  for  wave  propagation  in  single  materials. 


The  experimental  component  of  the  research  program  was  completed  as  envisaged 
in  the  proposal.  It  had  the  following  components; 

1.  Development  of  techniques  for  evaluation  of  material  parameters  that  appeared 
in  the  theoretical  models  considered. 

2.  Construction  of  a  uniform,  fully  saturated  sample. 

3.  Application  of  motion  to  the  sample  on  a  shaking  table. 

4.  Recording  of  input  acceleration  and  pore  pressures  up  to  liquefaction. 

5.  Analysis  of  the  experimental  data  for  the  purpose  of  evaluating  the  theoretical 
models. 

A  fine  to  medium  grained  sand,  Ottawa  sand,  was  chosen  as  the  material  to  be 
used  in  the  experimental  studies.  A  program  of  static  tests  aimed  at  identifying  basic 
material  parameters  that  appear  in  the  theoretical  models  was  performed.  For  liquefac¬ 
tion  experiments  particular  emphasis  was  placed  on  laboratory  techniques  for  proper 
preparation  of  the  necessary  sand  samples.  Suitable  methods  for  saturating  the  approxi¬ 
mately  SO  kilograms  of  sand  required  to  build  a  sample  were  developed  and  samples 
with  a  high  degree  of  uniformity  were  repeatedly  achieved.  Instrumentation  was  cali¬ 
brated  under  both  cyclic  and  static  loading  conditions.  Different  types  of  pore-pressure 
gage  applications  were  studied.  Initial  shaking  table  tests  were  designed  to  identify  the 
mast  reliable  and  sensitive  instrumentation  capable  of  recording  input  acceleration  and 
progressive  development  of  pore  pressures  up  to  liquefaction.  Improvements  in  collection 
of  data  and  its  processing  were  implemented.  Additional  instrumentation  was  acquired 
to  collect  an  increased  amount  of  information  during  tests. 
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These  tests  established  the  capability  of  the  experimental  set-up  to  serve  the 
immediate  needs  of  the  program.  Both  isotropically  and  anisotropically  consolidated 
nmplw  were  tested  to  liquefaction.  Input  motions  consisted  of  both  harmonic  accelera¬ 
tion  time  histories  and  random,  white  noise,  accelerations.  Further  tests  were  carried 
out  to  study  the  effect  of  overburden  an  a  soil  subjected  to  shaking.  The  experimen¬ 
tal  data  obtained  were  compared  with  the  results  of  previous  laboratory  investigations. 
Shortcomings  in  these  earlier  test  procedures  were  identified. 

Figure  2  illustrates  the  scope  of  the  actual  work  plan  and  the  accomplishments 
under  the  research  program. 

1.7  STRUCTURE  OF  THE  REPORT 

Section  D  contains  a  review  of  previous  theoretical  investigations  on  the  behavior 
of  saturated  soils  subjected  to  dynamic  loading  in  order  to  highlight  the  basic  assump¬ 
tions  of  the  various  approaches.  Section  ID  describes  a  dynamical  theory  of  motion  of 
saturated  soils,  developed  as  part  of  the  present  research  program,  based  on  description 
of  the  motion  of  a  given  set  of  soil  particles  in  convected  coordinates  and  regarding 
the  fluid  motion  as  relative  to  the  soil  particles.  Section  IV  describes  some  anlytical, 
and  semi-analytical  solution  procedures  developed  to  obtain  solutions  to  some  simple 
one-dimensional  problems  as  well  as  development  of  finite  element  techniques  for  anal¬ 
ysis  of  wave  propagation  in  saturated  soils.  Section  V  reviews  previous  laboratory 
investigations  and  summarizes  the  experimental  component  of  the  research  program. 
Section  VI  summarizes  the  results  of  the  research  and  contains  recommendations  for 
future  work.  Appendix  A  describes  some  concepts  related  to  the  use  of  convected  coor¬ 
dinates  to  describe  the  motion  of  a  solid.  Appendix  B  lists  the  research  publications  and 
presentations  that  resulted  from  the  research  effort. 


Fi&y«_2:  Research  Accomplished  Under  AFOSR  Grant  83-0055. 


Section  II 


REVIEW  OF  PREVIOUS  WORK 


2.1  INTRODUCTION 

Existing  theories  of  dynamics  of  fluid-saturated  soils  were  reviewed.  Limitations 
of  various  theories,  their  similarities  as  well  as  inconsistencies  were  identified.  The 
"engineering  approach”,  described  herein  and  based  on  methods  of  structural  mechanics, 
was  considered  along  with  the  continuum  mechanics  approaches.  Herein  we  briefly 
review  the  salient  features  of  these  approaches  as  the  background  for  the  theory  devel¬ 
oped  in  the  course  of  this  research.  Detailed  documentation  of  this  review  is  available 
in  the  technical  report  listed  as  item  1.2  in  Appendix  B. 

2.2  THE  ENGINEERING  APPROACH 

This  approach,  introduced  by  Seed  and  his  co-workers  [77-80,1481  uses  methods  of 
structural  dynamics  to  solve  the  problem  of  shear  wave  propagation  in  soils.  The 

approach  has  been  successfully  applied  to  several  case  histories  [85,88,128,145-154].  It 
consists  of  a  finite  element  analysis  of  the  dynamic  system  to  evaluate  the  stress  histo¬ 
ry.  This  is  followed  by  a  laboratory  study  of  the  material  behavior  under  cyclic  stress 
conditions  equivalent  to  those  determined  from  the  finite  element  analysis.  This 
approach  has  been  extended  to  include  a  periodic  updating  of  material  behavior  to 

allow  for  the  strain  history  as  well  as  pore- water  pressure  build-up  and  dissipation. 
The  stiffness  is  assumed  to  be  a  function  of  the  volumetric  strain  and  the  effective 

stress  in  the  soil.  Generation  of  pore  water  pressure  is  assumed  to  be  related  to  vol¬ 
ume  changes  in  water  and  soil.  Assuming  water  to  be  incompressible,  an  incremental 
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relationship  is  proposed  between  the  incremental  pore  pressure  and  the  volume  change 
in  the  soil.  For  one-dimensional  idealization  of  a  horizontally  layered  system  subjected 
to  shear  wave  at  the  base,  dissipation  of  the  pore  water  pressure  is  expected  to  be  gov¬ 
erned  by  the  equation 

'*Erl£,r..l..  +  EAk  (2-D 

where  w  is  the  unit  weight  of  water,  K  the  permeability,  E,  the  the  one-dimensional 
rebound  modulus  of  the  material  at  the  effective  stress  applicable  to  the  increment,  and 
etJ  are  components  of  the  strain  tensor  for  the  soil  Finn  [42]  called  this  the  coupled 
theory  of  liquefaction.  The  sequence  of  occurrences  is  assumed  to  be  as  follows:  shear¬ 
ing  stresses  cause  volume  changes,  volume  changes  result  in  pore-water  pressure  changes, 
pore-pressure  dissipation  follows,  pore-pressures  determine  effective  stresses  and  effective 
stresses  along  with  the  cumulative  shearing  strain  define  the  effective  shear  modulus  to 
be  used  for  determination  of  displacements  and  stresses  for  the  next  time  step. 

Item  1.4  in  Appendix  B  contains  details  of  the  methodology  of  the  "engineering 
approach”  to  liquefaction  as  well  as  its  implementation  in  a  computer  program.  The 
computer  code  developed  was  used  to  obtain  the  response  of  a  layered  soil  system  sub¬ 
jected  to  sinusoidal  base  acceleration.  The  problem  data  were  taken  from  a  case  study 
reported  by  Finn  [42].  The  surface  acceleration  and  displacement  as  well  as  the  time  to 
liquefaction  reported  by  Finn  and  obtained  using  the  code  were  in  good  agreement. 
However,  the  detailed  response  history  obtained  using  the  implicit  and  the  explicit 
methods  was  quite  different.  With  refinement  of  the  time  mesh,  one  would  expect 
the  explicit  method  to  yield  a  sequence  of  results  which  would  converge  to  the  solu¬ 
tion  obtained  by  the  implicit  method.  However,  the  solution  procedure  requires  experi¬ 
mental  data  which  are  dependent  upon  the  frequency  of  cycling,  the  natural  frequency 
of  the  soil  system  etc.  and  would  have  to  be  generated  for  each  case  studied.  The 
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study  showed  the  approach  to  be  cumbersome  and  based  upon  several  assumptions.  The 
physical  properties  required  as  data  for  the  method  can  only  be  determined  as  a  func¬ 
tion  of  the  complete  number  of  cycles  of  stress  at  a  certain  amplitude.  Thus,  in  an 
explicit  type  solution  scheme,  it  is  possible  to  reduce  the  time  interval  for  updates  to 
only  as  small  as  the  time  period  of  vibration.  For  this  reason  it  is  not  possible  to  gen¬ 
erate  a  convergent  sequence  of  solutions  based  upon  reducing  the  size  of  the  time  step. 
Implicit  schemes  are  expected  to  be  more  reliable.  Post-liquefaction  distribution  of  pore 
pressures  and  the  extent  of  the  resulting  displacements  in  the  system  cannot  be  cor¬ 
rectly  determined  in  this  theory.  Furthermore,  the  approach  has  only  been  used  for 
one-dimensional  wave  propagation,  and  cannot  be  readily  extended  to  two  and  three- 
dimensional  cases.  Use  of  this  theory  requires  considerable  experience  and  "judgement”, 
in  addition  to  extensive  laboratory  testing  program,  to  get  useful  results. 

2.3  CONTINUUM  MECHANICS  APPROACHES 

2-3.1  Biot's  Theory. 

Biots  theory,  based  on  concepts  of  coupled  motion  of  soil  and  water,  has  been  the 
most  popular  alternative  to  the  empirical  approaches.  This  pioneering  work,  was  based 
on  certain  postulates  regarding  description  of  motion,  notion  of  partial  stresses,  the  exis¬ 
tence  of  energy  and  dissipation  functions  for  the  saturated  mass  and  a  certain  form  for 
the  kinetic  energy  of  the  mixture.  These  assumptions  led  directly  to  the  conclusion 
that  constitutive  coupling,  inertial  coupling  and  equilibrium  couplings  exist  and  are 
symmetric  in  form.  The  equations  of  momentum  balance  were  written  for  the  mixture 
as  a  whole  and  for  fluid  motion  relative  to  the  solid.  Several  different  forms  of  Biot's 
theories  exist.  The  most  general  form  includes  body  forces,  inertia  forces  and  the  effects 
of  coupling  of  the  fluid  and  the  soil  mass  as  well  as  the  constitutive  coupling 


between  the  soil  and  the  water. 
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Biot  [ll]  wrote  the  constitutive  equations  for  the  flow  of  a  compressible  fluid 
through  a  porous  saturated  linearly  elastic  medium  assuming  the  existence  of  an  energy 
function  quadratic  in  measures  of  deformation  associated  with  the  saturated  soil  mass. 
These  included  the  six  components  of  the  soil  strain.  Change  in  water  content  was 
added  as  the  seventh  kinematic  variable.  While  extending  the  analysis  to  compressible 
fluids  and  anisotropic  elastic  or  viscoelastic  solids,  Biot  [14, IS]  introduced  the  volumetric 
strain  of  the  fluid  as  the  additional  strain  parameter  instead  of  the  change  in  water 
content  used  in  the  earlier  theory  and  explained  later  [20]  that  the  two  variables  were 
essentially  the  same.  Garg's  [50]  formulation  can  be  shown  to  correspond  to  Biot's.  In 
Garg's  theory,  the  constants  are  related  to  the  properties  of  the  constituents  and  their 
volume  fractions.  For  the  dissipative  case,  a  dissipation  function,  quadratic  in  relative 
velocity  was  introduced.  For  a  statistically  isotropic  saturated  material,  Biot  [17-19] 
expected  the  kinetic  energy  to  be  quadratic  in  the  velocities  of  the  fluid  and  the  soil 
and  a  coupling  term  was  included.  This  introduced  an  inertial  coupling  between  the 
soil  and  the  fluid.  However,  it  is  difficult  to  assign  numerical  values  to  the  various 
quantities  that  arise  as  a  result  of  this  coupling.  As  a  part  of  this  review,  the  theory 
was  implemented  in  a  finite  element  computer  program  and  a  parametric  study  carried 
out  to  investigate  the  efect  of  this  coupling  on  the  dynamic  response.  Preliminary  stud¬ 
ies  indicated  that  the  effect  would  be  insignificant.  However,  further  study  of  this 
particular  feature  is  needed. 

While  developing  finite  element  solution  procedures  for  the  problem,  Ghaboussi 
[53],  following  Biot  [15,18,191  introduced  relative  volumetric  strain  in  the  formulation. 
The  momentum  balance  and  the  continuity  equations  were  written  in  terms  of  six  dis¬ 
placement  components  viz.  the  soil  displacements  and  the  relative  displacements  of  the 
fluid.  A  Rayleigh  type  viscous  damping  term  was  introduced.  This  intrinsic  clamping  of 
the  soil  is  in  addition  to  the  damping  associated  with  relative  motion.  Ghaboussi  [53] 
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developed  a  variational  formulation  for  the  problem  but  for  the  purpose  of  finite  ele¬ 
ment  analysis  he  used  the  Galerkin  procedure.  Also  the  boundary  conditions  were  not 
treated  properly. 

133  Theories  of  Mixtures. 

23.2.1  Introduction 

To  apply  the  principles  of  continuum  mechanics,  it  is  customary  to  regard  a 
fluid-saturated  solid  as  superposed  continua.  Averaging  of  various  quantities,  kinematic 
as  well  as  mechanical,  associated  with  the  various  constituents,  is  inherent  in  this 
assumption.  The  mixture  is  defined  by  the  current  coincident  configuration  of  the  con¬ 
stituents.  It  is  assumed  that,  in  the  current  configuration,  each  point  of  space  is  occu¬ 
pied  by  a  particle  of  each  of  the  constituents.  This  necessitates  the  introduction  of 
"bulk”  description  of  the  material  instead  of  the  "intrinsic”  description  which  would 
apply  if  the  material  were  the  single  constituent  of  a  body.  In  any  theory  of  mix¬ 
tures,  it  would  be  necessary  that  as  the  volume  fraction  of  one  of  the  constituents 
approaches  unity,  and  the  remaining  constituents  disappear,  the  theory  for  a  single  con¬ 
stituent  be  realized  as  a  limiting  case. 

In  developing  a  rational  theory  of  mixtures,  Truesdell  [165]  laid  down  the  fol¬ 
lowing  principles: 

1.  All  properties  of  the  mixture  must  be  mathematical  consequences  of  properties 
of  the  constituents. 

2.  So  as  to  describe  the  motion  of  a  constituent,  we  may  in  imagination  isolate 
it  from  the  rest  of  the  mixture,  provided  we  allow  properly  for  the  actions 
of  the  other  constituents  upon  it. 

3.  The  motion  of  the  mixture  is  governed  by  the  same  equations  as  is  a  single 
body. 

The  third  of  these  "principles”  is  open  to  serious  objection.  This  issue  is  addressed  in 
a  later  paragraph. 
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2322  Density  of  Each  Constituent  and  of  the  Mixture. 

The  effective  (also  referred  to  as  partial  or  bulk)  density  of  the  kth  constituent, 
when  it  is  regarded  as  one  of  the  continua  occupying  every  point  in  the  spatial  region 
of  interest,  is  related  to  its  intrinsic  density  as 

(22) 

where  nw  is  the  volume  fraction,  superscipted  Roman  characters  enclosed  in  parentheses 
identify  the  constituent  in  the  mixture,  and  a  superscripted  asterisk  denotes  an  intrinsic 
quantity.  The  density  of  a  mixture  of  n  constituents  is  defined  as  the  weighted  sum 
of  the  densities  of  the  constituents  ien 
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P  =  2>tt’  (2.3) 

k-1 

Other  definitions  for  partial  densities  have  been  used.  Terzaghi  [161]  assumed  the 
effective  density  of  water  to  have  the  intrinsic  value  and  regarded  the  soil  as  buoyant 
in  the  water. 


23.23  Description  of  Motion. 

Several  approaches  have  been  used  to  describe  the  motion  of  the  constituents  of  a 
multicomponent  mixture.  Often  the  deformation  is  referred  to  an  initial  configuration 
for  each  constituent  and  motion  to  the  place  coordinates.  The  equations  of  balance  are 
written  for  a  fixed  volume  in  space.  Another  approach  is  to  refer  the  motion  of  all 
constituents  to  the  reference  configuration  of  one  of  them.  Yet  another  is  to  refer  all 
motion  to  the  current  configuration  which  is  the  same  for  all  constituents.  Superposi¬ 
tion  of  relative  diffusive  motion  of  the  constituents  upon  the  mean  motion  of  the 
mixture  as  a  whole  is  also  used.  For  a  binary  mixture  of  a  solid  and  a  fluid,  some 
investigators  describe  the  motion  of  the  solid  with  respect  to  its  reference  configuration 
but  the  motion  of  the  fluid  is  described  as  relative  to  the  solid.  Another  approach  is 
to  refer  to  a  material  region  consisting  of  the  same  set  of  particles  of  one  of  the  con- 
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stituents  so  that  the  bounding  surface  of  this  constituent  varies  with  time.  If  the 
description  of  the  material  particles  is  the  same  during  motion  as  in  the  initial  config¬ 
uration,  this  would  be  the  convected  coordinate  description  for  the  reference  constituent. 
Gibson's  theory  of  consolidation  [SS]  is  a  special  case  of  this  approach. 


Using  a  fixed  rectangular  cartesian  reference  frame,  the  deformation  gradient  is 
defined  as  the  partial  derivative  of  the  place  coordinates  with  respect  to  the  reference 
coordinates.  The  velocity  vectors  and  the  acceleration  vectors  are  defined  as  the  partial 
time  derivatives  of  the  place  coordinates  for  the  same  material  particle.  Most  investiga¬ 
tors,  also  introduce  a  barycentric  velocity  rate  for  the  mixture  as  a  whole  as: 


v. 


_  n  p(k)v(» 

i  Z*  n 


(2.4) 


k-1 


where  v,  are  components  of  the  velocity  vector.  A  material  rate  for  the  mixture  is 
established  using  the  identity: 


Dt 

k-l 


p-5. 

Dt 


(2.5) 


The  diffusive  velocity  is  defined  as  the  velocity  of  a  constituent  relative  to  that  of 
the  mixture,  i-e, 

uk)  =  v[k)-  v.  (2.6) 

Biot  [15,18,19],  for  the  case  of  binary  mixture  of  a  solid  and  a  fluid,  introduced  a 
nominal  relative  velocity 

w,-  (2.7) 

This  and  its  variants  were  used  by  Ghaboussi  [53],  Krause  [86]  and  Kenyon  [84],  among 


others. 


18 


23.2^4  Measures  of  Deformation. 

The  deformation  gradient  and  the  spatial  velocity  gradient  can  be  used  to  com¬ 
pletely  define  rates  of  deformation.  For  material  objectivity  to  be  satisfied,  the  vorticity 
tensors,  in  binary  mixtures,  must  occur  as  the  difference  of  the  vorticities  of  the  con¬ 
stituents. 


For  a  fluid-saturated  porous  solid,  the  deformation  gradient  can  be  used  to 
describe  the  changing  configuration  of  the  solid.  Garg  [48,491  and  Morland  [103,104] 
among  others,  used  as  measures  of  deformation  of  the  solid,  the  quantity 

a*) 

where  u,  denotes  components  of  the  displacement  vector.  Here,  and  in  the  sequel,  we 
use  the  standard  indicial  notation.  The  Roman  indices  take  on  values  in  the  range  1, 
2,  3.  Summation  on  repeated  indices  is  understood  unless  stated  otherwise.  A  subscript¬ 
ed  comma  denotes  differentiation  with  respect  to  the  coordinate  represented  by  the  sub¬ 
script  following  it.  Parentheses  around  a  pair  of  indices  indicate  "symmetric  part"  and 
square  brackets  the  "antisymmetric  part”.  Superposed  dots  indicate  differentiation  with 
respect  to  the  time  parameter.  This  measure  of  deformation  characterizes  the  linear 
theory.  For  setting  up  constitutive  relations  for  flow  through  a  porous  saturated  elastic 
anisotropic  medium,  Biot  [15,16]  used  relative  volumetric  strain  defined  by  the  equation 

$  =  n(2)(e(k2>-e£)  (2.9) 

Other,  more  general  measures  of  deformation  are  also  used,  e.g.  [103].  For  an  initially 


isotropic  solid,  it  is  possible  to  define  as  a  measure  of  deformation  the  quantity 


2e.  =  F  F  -8 

ij  mi  mj  ij 


(2.10) 


where  Fu  are  components  of  the  deformation  gradient  tensor.  If  the  only  deformation 
of  interest  is  the  volume  change  the  change  in  density  is  an  adequate  measure  of 
deformation.  In  a  binary  fluid-solid  mixture,  for  small  strain 
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„(D  rt(D _ -d)  (i) 

P  -Pa - P  eH 


(2.11) 


where  the  subscript  denotes  initial  value  or  the  value  in  the  reference  configuration. 
For  the  fluid  constituent  the  deformation  is  completely  represented  by  its  density,  uen 


(2)  (2)  (2Y 

(2)  Pq  _  1  _  no  Pq  _  | 
»  ~  n‘V21- 


Introducing  an  intrinsic  measure  of  deformation  [103] 

Morland  [104]  also  proposed,  for  infinitesimal  strain  theory, 
Cji  "  [l-n<2)] 


(2.12) 


(2.13) 


(2.14) 


— e(i>*  s  dl  -i 
P(o1} 


(2.15) 


e^-ie<»,8..-e!.,>-ie(,)  8. 

ij  3  mm  ij  ij  J  mm  ij 


(2.16) 


and  [106] 


e(2)"  =  e(2)  + 
jj  jj 


t  n!2>- 


(2.17) 


In  addition  to  the  measures  of  deformation  defined  above,  the  volume  fractions  of  the 
constituents  have  been  used  as  state  variables.  An  intrinsic  rate  of  strain  of  the  fluid 


d<2>-  _  d<2>+  I  ^o_  U 


Dt  n<«  'j 


(2.18) 


was  introduced  by  Morland  [104].  Here  —  denotes  the  "material  rate”.  Another 

measure  of  relative  deformation  sometimes  used  is,  (e.g,  Krause  [86D,  the  change  in  flu¬ 
id  content  of  a  fixed  volume  in  space.  This  amounts  to 

_  (2)  (2)  (  (2)  (2K  (2f 

n  =  p  -p0=(u  - n 0  )p 


(2.19) 


Carroll  [27],  considering  a  fixed  volume  of  the  solid  in  the  reference  configura¬ 
tion,  assumed  the  total  deformation  to  be  the  sum  of  the  strains  of  the  solid  particles 


and  the  change  in  geometry  of  the  pores.  The  bulk  volume  strain  was  found  to  be 


(2) 


AY  Anw  AV<1) 
n(I) 


V  n<U  ' 

For  the  anisotropic  case,  Caroll  [29]  wrote 


(2.20) 


•b-'T+- 


(2.21) 


as  the  strain  of  the  reference  material  volume  of  the  bulk  solid.  Carroll  [28]  assumed 
the  solid  intrinsic  deformations  to  be  reversible  under  solid  pressure  and  the  pore  vol¬ 
ume  change  to  be  irreversible  under  effective  pressure.  This  would  explain  cumulative 
volume  change  under  cycling  of  load.  Kenyon  [84]  regarded  the  specific  volume  of  the 
solid  as  a  state  variable  related  to  the  intrinsic  pore  fluid  pressure  and  the  solid  defor¬ 
mation  gradient. 

23.2~5  Balance  Laws 


a).  Earlier  Theories. 


The  equation  of  mass  continuity  of  the  fluid,  equating  the  inflow  into  an  elemen¬ 
tary  volume  of  a  rigid  porous  one-dimensional  solid  with  the  increase  in  fluid  content, 
is 

4-(  P(2)v  )  +  £(  P<2> 4t-)  =  0  (2-22) 

qA  gt  gX 

Here  v  is  the  velocity  of  the  fluid  relative  to  the  soil  and  x,  X  are  the  coordinates  in 
the  current  and  the  reference  configurations,  respectively,  referred  to  a  Cartesian  refer¬ 
ence  frame.  Allowing  for  compressibility  of  the  pore  fluid,  Gibson  [55]  used  Scheideg- 
ger's  [140]  formulation  of  d’Arcy's  law.  The  flow  equation  has  the  form 
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Ghaboussi  [54]  stated  Biot's  momentum  balance  equations  as 
tijiJ+pbi  =  piij»-p®[  uS^-u^  =  p<l)u;i)+p(2)u<2) 


(224) 


and 


^,,+A  =  02s 

■where  k  is  a  measure  of  permeability  of  the  soil  and  ty  are  components  of  the  "total” 
stress  tensor.  Garg  [50]  wrote  Biot's  equations  of  motion,  in  the  absence  of  body  forces, 
in  the  following  form: 

^1-p<,)S<,>+b(i;,>-i«a)  (2.26) 

IT,  =  (2.27) 

The  second  term  on  the  right  hand  side  of  these  equations  represents  the  viscous 

coupling  between  the  solid  and  the  fluid. 


b).  Mechanics  of  Mixtures. 

For  each  constituent  of  the  mixture,  Truesdell  [164-166]  postulated  the  balance  laws 
in  the  point  form.  Using  the  notation  of  [64],  these  were: 

i.  Continuity  of  Mass. 


(2.28) 


where  cU)  is  the  mass  production  fraction  of  the  kth  constituent.  An  alternative  form 
of  the  mass  continuity  equation  is 


D(k)  „<k)j.  «(kUk>_  ~.(k) 

— p  +p  V,  =pc 


Bowen  [24]  wrote  the  above  equation  in  the  form 
^Up(k,detFj>)  =  pc(k)detlFj)l 


(2.29) 


(2.30) 


I 
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iL  Balance  of  Linear  Momentum. 


Acp^vi11)  ♦  ( A!' v;1),- ,« = *»;■’+  p(Ub“ 


(2.31) 


Here  t^.pb^.pm^  are,  respectively,  components  of  the  partial  Cauchy  stress  tensor,  the 
body  force  vector  per  unit  mass  and  the  partial  momentum  supply  density  for  the 


constituent.  An  alternative  form  is: 


p<k)f*k)  _  t<k)j+  pm|k)+  p(kVk)—  p  c(k Vk 


(2.32) 


Bowen  [24]  proposed  to  replace  the  momentum  supply  term  m^  by  m,(k)+  c^v^  to  get 


p(k¥k)  =  t<k)  +  pm<k)+p(kVk) 


(2.33) 


iiL  Balance  of  Moment  of  Momentum. 


(234) 


Here  pM^  are  components  of  the  skew-symmetric  tensor  describing  the  partial  produc¬ 
tion  density  of  the  moment  of  momentum  of  the  kth  constituent. 


Iv.  Balance  of  Energy  Rates. 


p««  D_(  u«)  +  q“-  +  »4U'U 


(2.35) 


where 


Al/k)  =  e(k)— ( l/k>-  Iv(kVk>)  c(k)-  v'k)m(k>-  IM<k)v(k) 

2  1  1  11  2  '•> 


(2.36) 


Here  l/u ,  qU),  r<4),  e(k)  are,  respectively,  the  specific  internal  energy,  components  of  the 
partial  heat  flux  vector,  the  partial  energy  supply  and  the  partial  energy  production 
density  of  the  kth  constituent 
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Bowen  [24]  wrote  the  energy  rate  balance  equation,  for  a  fixed  volume  in  space 
of  the  k.th  constituent,  and  derived  the  point  form  of  the  equation  of  energy  balance 

_<k)D(k)l/k)  ,  _(k)  Jk)_(k)_  .(k)„(k),  „0 

p  -DT”+qi -i  P  f  _tMV‘.J+pe 


(k) 


(2.37) 


Several  alternative  forms  of  the  above  equation  have  been  stated  [24,164].  Bowen  [24] 
pointed  out  that  for  the  case  of  single  temperature  mixtures,  explicit  use  of  an  energy 
equation  for  the  constituents  is  not  needed. 


Truesdell  [164,166]  postulated  balance  equations  for  the  mixture  as  well.  These 
equations  can  be  derived  by  summing  over  the  equations  for  the  constituents.  Truesdell 
introduced  the  following  quantities  to  ensure  that  the  resulting  equations  for  the 
motion  of  the  mixture  have  the  same  form  as  the  equations  of  motion  of  a  single  con¬ 
stituent. 


L  Specific  internal  energy  of  the  mixture 

p  **  2  '  * 

r  k- 1 

Or,  equivalently, 

p(  u + iv.v.)  =  £P(k)(  u<k) + iv;kVk)) 

k-1 


ii.  Total  stress  tensor 

t  =  f(t(k)-p(k)u(k)u<k>) 

ij  L*  ij  r  i  j 

k-l 

Or,  equivalently, 

t  -pvv  =  T(  tjk)-p(k)v(k)v(k>) 

ij  r  I  J  t-t  ij  r  I  J 


k-l 


(238) 


(239) 


(2.40) 


(2.41) 
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iii.  Total  heat  flux  vector 


q,  -  £<*"+*<**  U"+  luf 


k-I 


Or,  equivalently. 


(2.42) 


q,+  pv,[U+ivJvJ]-tljvJ 


=  £[*"+ Af'tu"* 

k-I 


(2.43) 


iv.  Specific  energy  supply. 

r  =l£p<«(  r‘u+b!uujl>)  (2.44) 

**  k-1 

Or,  equivalently, 

p  (  r  +  b.  v.)  =  £  p(k)(  r(k)+  bjkVk>)  (2.45) 

k-l 

Also,  thermodynamic  isolation  of  the  mixture  was  assumed,  Le, 


Q 


o 

tl 

j5 

'~b 

(2.46) 

k-l 

ixk) - o 

(2.47) 

k-l 

£yk)  =  o 

ij 

k-l 

(2.48) 

£e(k)  =  0  (2.49) 

k-l 

Kelly  [82i  Truesdell  [164],  Green  [66]  and  Bowen  [24]  wrote  the  equations  of 
mass,  momentum  and  energy  conservation  for  the  mixture  contained  in  a  volume  V 
bounded  by  an  arbitrary  fixed  surface  A.  The  conservation  of  mass  was  expressed  by 


JL 

dt 


dV  + 


/iv 


%“dA 


A  k-l 


=  0 


(2.50) 


25 


Equating  the  linear  change  of  momentum  to  the  total  force  exerted  on  the 
material,  using  the  divergence  theorem.  Green  [66]  obtained 


|  p£v +||/W 


(2.51) 


The  point  form  immediately  follows  from  this  integral  form  of  the  equation  of  bal¬ 
ance  of  momentum.  Green  [66]  noted  that  with  the  definition  of  the  total  stress  tensor 
proposed  by  Truesdell,  their  form  of  the  equation  is  the  same  as  that  obtained  by 
Truesdell.  However,  if  the  total  stress  is  defined  as  the  sum  of  the  partial  stresses,  the 
equation  reflects  the  fact  that  the  total  rate  of  increase  of  linear  momentum  is  not 
equal  to  the  barycentric  rate  of  increase  of  momentum  of  a  continuum  of  density 
moving  with  the  barycentric  velocity. 


The  theory  for  a  mixture  of  two  constituents  presented  by  Green  [62,64]  was 
generalized  by  Mills  [99]  to  the  case  of  multicomponent  mixtures.  Green  [62,64]  consid¬ 
ered  the  concepts  of  stress,  heat  flux,  and  energy  supply,  assumed  to  be  primitive  to 
each  constituent,  to  be  primitive  to  the  mixture  as  a  whole  as  well.  It  was  proposed 
that  the  total  stress  and  the  total  heat  flux  for  the  mixture  should  equal  the  sum  of 
the  corresponding  quantities  for  the  constituents.  Green  [66]  stated  that  TruesdelTs  equa¬ 
tions  were  correct  but  there  was  difficulty  accepting  the  interpretations  associated  with 
some  of  the  quantities  which  occur  in  these  equations.  For  a  volume  V  enclosed  by  a 
fixed  surface  A,  following  TruesdelTs  contention  that  the  rate  of  energy  equality  has 
the  same  form  as  for  a  single  constituent.  Green  [66]  wrote 


A 

6t 


p(  U  +  -jVv.)dV  +  J pnv.(  U  +  yVjV.JdA 
=  f£pk( r(k)+b'k)v[k))dV+  [(  Vj— q,)n,dA 

*V  k-1  JA 


(2-52) 
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Here  the  left  hand  side  represents  the  rate  of  change  in  energy  in  volume  V  bounded 
by  a  fixed  surface  A  plus  the  energy  flux  for  the  mixture  across  the  boundary.  U  is 
the  specific  internal  energy  per  unit  mass  of  the  mixture  and  is  related  to  internal 
energies  of  the  constituents  by  Truesdell's  equation.  The  equality  may  be  regarded  as 
written  for  a  surface  moving  with  a  continuum  which  has  a  velocity  field  equal  to 
the  barycentric  velocity.  Then  the  left  hand  side  is  equal  to  the  material  rate,  executed 
on  the  mixture,  of  the  sum  of  the  internal  energy  and  the  kinetic  energy  of  the  mix¬ 
ture.  However,  this  line  of  thought  is  open  to  objection.  It  would  assume  the  existence 
of  mixture  particles  and  the  time  rates  executed  on  them.  This  is,  in  general,  not  cor¬ 
rect  as  the  barycentric  velocity  is  not  particle  velocity  in  the  ordinary  sense.  Green 
[66]  accepted  the  form  of  (2-52)  but  not  the  interpretation  associated  with  some  of  the 
quantities  occuring  in  it. 


Green  [62]  proposed  a  rate  of  energy  equality  in  the  following  form; 

/<  u  £<■'■  v;1 + i  £p,iiv;iiv;iv»)  nj « 

®  *v  k-1  *A  k-1  k-1 


(  p  r  +  £p(k)bJkVk>)  dV  +  f£(  t'k Vk)-  q  )  dA  (2.53) 

k-1  *X  k-1 

Here,  the  heat  fluxes  and  the  energy  supply  are  assumed  to  be  additive.  U  is  the 
internal  energy  per  unit  mass  of  the  mixture  allowing  for  all  interactions  between  the 
constituents.  The  t,u)  are  components  of  tractions  associated  with  the  kth  constituent 
and  the  surface  A.  In  this  theory,  multipolar  stresses  and  externally  supplied  multipo¬ 
lar  body  forces  were  excluded.  Green  [62]  also  made  no  attempt  to  define  the  internal 
energy  and  the  energy  supply  for  each  constituent.  It  was  considered  unnecessary  for  a 
complete  general  theory. 
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Later,  [64l  the  role  of  interactions  between  constituents  was  made  explicit  by  writ¬ 
ing,  for  each  constituent,  the  rate  of  energy  equality 


(2.54) 


where  are  the  internal  force  and  couple  acting  on  kth  constituent  due  to 

interactions  and  q00  represent  volume  and  surface  contributions,  respectively,  to  the 
balance  of  energy.  Also 


n 


o 

II 

' 3 

w«f 

(2.55) 

k-1 

Q 

o 

II 

A  =• 

(2.56) 

k-i 

n 

(237) 

k-1 

Without  loss  of  generality,  AtJ  can  be  taken  to  be  [64]  antisymmetric  as  is  anti¬ 
symmetric.  Bowen  [24],  and  Bedford  [9]  used  essentially  the  same  form  in  writing  ener¬ 
gy  balance  equality  for  each  constituent.  Assuming  the  quantities  ,  p,  U,  v,  for  the 
mixture  are  sufficiently  smooth  in  the  space  and  time  variables,  it  can  be  shown  that: 

[p~-  -  pr  +  i>(kVk)(  ff- b<k>)]  dV  +  f  £pc(k,(U  +  iv<kVk^V 


=  f£(  t<kVk>-  q  )  dA  (2J58) 

*A  k-1 

Invariance  of  the  energy  equality  under  superposed  uniform  translational  velocities 
yields,  for  arbitrary  V, 
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1 

£  £pc(k)dV  =  0  (2-59) 

ia,  the  mass  elements  of  the  mixture  are  conserved.  It  was  also  shown  that 

• 

(2.60) 

1 

A  generalization  of  Cauchy's  stress  principle  was  derived  in  the  form 

-a 

if-  £$s 

k-l  k-l 

(2.61) 

In  point  form,  the  equation  of  balance  of  linear  momentum  for  the  mixture  postulated 

• 

by  Truesdell  was  obtained.  Green  [62]  obtained,  for  a  binary  mixture. 

• 

+  y  J* (tj°-  tj2>) (v<°-  vj2>)  dA  =  J  q  dA 

(2.62) 

Defining,  following  Mills  [99] 

Jkl  -  t<k)-r(k)n 

Pi  ~  ^  nj 

(2.63) 

• 

and  applying  the  rate  of  energy  equality  to  a  tetrahedron  bounded  by  the  coordinate 

planes  and  a  plane  with  unit  normal  nj  for  heat  flow  hj  across  plane  Xj,  for  a  binary 

• 

mixture.  Green  [62]  obtained: 

(2.64) 

In  the  special  case  when  q  =  hjn?  p,(1)  and  p^  vanish  except,  possibly,  in  the  case  of  no 

• 

* 

relative  motion  between  the  constituents.  Rate  of  energy  equality  can  also  be  stated  as 

f\p^r-pr+Kr  i*Mk>+ I KW  -  £<Mk: 

*V  l  k-l  k-l  k-l 

dV  =  0  (2.65) 

• 

where 

0<k)  = 

(2.66) 

• 

• 
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Invariance  under  a  superposed  uniform  rigid  body  angular  velocity  leads  to  symmetry 
of  the  sum  of  partial  stresses.  The  partial  stresses  do  not  have  to  be  symmetric.  The 
rate  of  energy  equality  was  written  in  point  form  as: 


DU 

Dt 


-  o  CM» 


k-1 


k-l 


k-1 


where 


p^f-iav^-v^ 


=  + p(k)c  b<k)-  ^  -  jc(k)c  v;k)-  v!d>) 


(2.68) 

(2.69) 


In  setting  up  the  rate  of  energy  equality,  it  was  stated  that  the  internal  energy 
per  unit  mass  of  the  mixture  may  not  be  equal  to  the  sum  of  internal  energies  of  the 
constituents.  Green  [66]  showed  that  if  the  sum  of  internal  energies  of  the  constituents 
is  defined  by 


pU*  =  J^p(k)U(k)  (2.70) 

k-1 

the  rate  of  energy  equality  in  the  form  of  equation  (2.58)  is  realized  if  the  energy  of 
the  mixture  is  defined  by  the  equation 


P 


DU 

Dt 


(2.71) 


where 


K  =  £(p(k)u(ik)U(k>)i  (2.72) 

k-1 


c).  Other  theories. 


Several  investigators  have  used  the  basic  concepts  introduced  by  Truesdell  and  by 
Green.  These  efforts  aim  at  simplifying  the  description  of  motion  for  certain  special 
cases.  For  instance,  in  granular  porous  media,  the  total  deformation  can  be  viewed  as 
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made  up  of  two  parts;  one  related  to  deformation  of  the  solid  particles  and  the  other 
to  their  rearrangement  Le.  changes  in  pore  geometry.  For  compressible  materials,  volume 
fractions  have  been  introduced  as  additional  variables.  Herein  we  outline  some  of  the 
results. 


L  Mass  Continuity  Equation  Using  Relative  Velocity 


n(k)+  (  n(kV.k))  j  =  0 


Krause  [86],  referring  to  fixed  volumes  in  space,  for  continued  saturation  and  no 
mass  production,  wrote  the  mass  continuity  equation  as 

A‘W  -  0  (2.73) 

If  the  materials  are  intrinsically  incompressible,  p^*  =  0.  Hence,  if  the  intrinsic  density 
is  also  spatially  constant, 

(2.74) 

For  a  binary  mixture,  the  above  equations  lead  to 

[  v(.°+  n(2)(  v(.2)-  v\l>)  ]  l  =  0  (2.75) 

Hsieh  [76]  added  the  equations  of  continuity  of  mass  for  each  of  the  two  constituents 
in  a  binary  mixture  and,  for  no  chemical  reaction,  obtained 

p  +  (  pvj °) I .  +  (  p(2) w.)  ,  =  0  (2.76) 

as  the  mass  continuity  equation  in  terms  of  relative  velocity.  Assuming  small  defor¬ 
mations,  he  also  wrote  the  continuity  equation  for  the  fluid  volume  contained  in  a 
fixed  volume  in  space  as 

JL  ( p(I)- p<0'>)  +  p^V  (p<02)  J_(u'2)-  u^)  =  0  (2.77) 

Integrating  over  time,  for  compressible  fluid  of  spatially  uniform  initial  density 


P  =  Pod-e^  +  Po’  J  w.  .Wdr 


(2.78) 
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ii.  Mass  Balance  in  Terms  of  Porosity 


Fukuo  [46]  used  the  equations  of  mass  balance  to  set  up  equations  in  terms  of 
volume  fractions  of  the  constituents.  Assuming  no  chemical  interaction,  using  intrinsic 


densities,  he  obtained: 


*(k)  (k)* .  («•(«• .  /  (k)  (k>  (kk  _  A 

np  +  n  p  +  ln  p  v.  1=0 
If  the  kth  constituent  is  incompressible,  p°°*  =  0,  and  p^*  =  0.  Hence, 
n(k>+  (  n(kVk>) .  =  0 


(2.79) 


(2.80) 


Hsieh  [76]  considered  a  porous  solid  saturated  with  an  incompressible  fluid  and  under¬ 
going  small  deformations.  For  this  condition,  considering  a  unit  volume  in  the  unde¬ 


formed  state,  they  showed  that  the  fluid  content  change  is 


n  =  -(n^-n^-n^1 


They  also  showed  that 

$t  °  °  a t  ij 


(2.81) 


(2.82) 


This  is  a  relationship  between  rate  of  porosity  change,  the  rate  of  volumetric  strain 
and  the  relative  velocity  vector.  These  quantities,  in  a  theory  for  incompressible  fluids 
and  no  thermal  or  chemical  effects  cannot,  therefore,  be  treated  as  independently  vari¬ 
able.  For  compressible  fluids,  Morland  [103]  proposed  constitutive  equations  for  porosity. 


iii.  Alternative  Form  of  Linear  Momentum  Balance. 


Considering  balance  of  momentum  of  fixed  volumes  in  space,  Hsieh  [76]  derived 
the  local  form  of  the  momentum  balance  equation,  for  no  body  forces,  as: 


f  P,k)v(k)+  f  ( p(k)v(kVk>)  =  y  t(k) 

At  “  ‘  J  ■]  £*}>.} 


(2.83) 


Assuming  additivity  of  stresses,  for  a  binary  mixture,  he  obtained 
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p  -5-  v(1)+  p(2^  -2-  v  +  v(2)v  +  v  v  J  =  t  (2  84) 

PDt  *  P  L  1  J  ».J  i  i.r  j».j  K  J 

where  v,  =  —  v[°  is  the  relative  velocity. 

iv.  Energy  Balance  in  Terms  of  Porosity. 

Goodman  [57]  postulated  the  equation  of  energy  balance  for  a  porous  material 
with  porosity  n  as 

/*np‘|(U  +  -iv.v.  +  IKT45-)  )~b.v  — il°--r  dV 
Dt  4  P[  2  “  2  ^3t^  *  *  at 

=  J  (tjiv.  +  sil2.-qi)n.dA  (2.85) 


Here  K,  s,  1  are,  respectively,  the  equilibrated  inertia,  components  of  the  equilibrated 
stress  vector,  and  the  external  equilibrated  body  force.  This  equation  admits  an  addi¬ 
tional  degree  of  freedom,  viz,  the  volume  fraction.  A  kinetic  energy  term  was  associat¬ 
ed  with  the  rate  of  change  of  n.  Similarly,  rate  of  work  terms  were  associated  with 
the  rate  of  change  of  n  over  volumes  and  surfaces  using  generalized  forces  Sj  and  L 
Invariance  of  this  equality,  as  in  Green  theory,  leads  to  the  equations  of  linear  momen¬ 
tum  and  angular  momentum  balance.  Goodman  [57]  also  postulated  an  equation  of 


motion,  called  the  equation  of  balance  of  equilibrated  forces,  for  the  variable  n,  as 


£{“pUaT+1+gldV".(w1A 


(2.86) 


Here  g  is  the  intrinsic  equilibrated  body  force.  They  also  wrote  the  local  form  of  this 
equation. 

v.  Other  Form  of  Energy  Balance  Equation  Bowen  [24]  postulated  the  point 
form  of  the  rate  of  energy  equality  as 


p  [  U  +  —  v  v  ]  =  ( t  v  —  q.)  +  p  r  +  7p(k,v(l> b(k) 
2  '  '  J»  J  I*  •  ' 


(2.87) 
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Using  an  interaction  force  vector.  Green  [67]  wrote  the  energy  balance  equation  in 


point  form  as 


(2.88) 


where 


*  =  £( 


(2.89) 


pu  =  £P(k)u(k) 


(2.90) 


=  t^+Ab^-^-pc^ 


(2.91) 


Mokadam  [lOl]  considered  the  total  internal  energy  for  the  mixture  to  consist  of 


three  components. 


E  =  U  +  T  +  L 


(2.92) 


where  U,  T,  L,  are  respectively,  the  molecular,  the  kinetic,  and  the  potential  energies 
per  unit  mass  of  the  mixture.  Identifying  the  diffusive  force  as  the  body  force  causing 
mass  flow,  Mokadam  postulated  equations  of  balance  of  energy  in  the  form 


*  -(,,Uvi+l'r,iivi).j+Djv) 


(2.93) 


where  Dj  are  components  of  the  diffusive  force.  Also 


(2.94) 


For  no  chemical  reaction,  the  above  equation  along  with  mass  conservation  implies 


£L=fv 

Dt  1  1 


(2.95) 


Mokadam  wrote,  for  no  chemical  reaction 


A ^  =  ( t  v  )  — ( pTv  +  q  ) 

a  t  <j  i  .  j  r  i  ,i 


(2.96) 
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Constitutive  Relations 

In  order  to  set  up  constitutive  equations,  it  is  necessary  first  to  define  the 

mechanical  quantities  for  which  such  relationships  are  desired  and  to  identify  the  kine¬ 
matic  or  state  variables  on  which  these  quantities  might  depend.  For  a  mixture  of  sev¬ 
eral  constituents,  the  stresses  in  the  constituents  are  obviously  the  primary  mechanical 
variables.  Truesdell  [163]  following  Maxwell  [97],  recognized  diffusive  resistance  as 

another  mechanical  variable  reflecting  the  interaction  between  the  constituents.  The  rate 
of  energy  equality  contains  scalar  products  of  "corresponding”  quantities.  This  indicates 
the  need  for  postulating  constitutive  relations  for  components  of  the  symmetric  and  the 
antisymmetric  parts  of  the  partial  stress  tensors  and  an  interaction  quantity.  Again,  if 
the  equations  of  mass,  linear  momentum,  angular  momentum,  and  energy  balance  for 
constituents  be  regarded  as  equations  for  certain  quantities,  constitutive  equations  are 

required  for  the  other  quantities  appearing  in  those  equations  and  also  for  the  partial 
entropy  of  each  constituent.  These  constitutive  equations  are  subject  to  the  balance 
equations  for  the  mixture  and  to  an  appropriate  entropy  production  inequality. 

There  has  been  some  difficulty  in  defining  components  of  the  partial  stress  ten¬ 
sors.  For  fluid-saturated  solids,  the  isotropic  fluid  stress  is  generally  considered  to  be 
the  stress  variable  in  addition  to  the  stresses  in  the  solid.  Biot  [ll]  regarded  the  total 
stress  and  the  fluid  stress  as  the  mechanical  variables.  The  definition  of  pore-fluid 

pressure  used  by  various  investigators  differs  considerably.  Traditionally,  for  a  water- 
saturated  soil,  the  pressure  recorded  by  piezometers  inserted  into  the  water-filled  pores 
has  been  assumed  to  be  the  fluid  pressure  acting  over  100  percent  of  the  area  of  inter¬ 
nal  surfaces  [5,26,74,81,90,162,175].  Biot  [17]  pointed  out  that  the  generalized  forces 
defined  by  divergence  of  the  stresses  are  correctly  defined  by  the  virtual  work  of 
microscopic  stresses  per  unit  value  of  the  displacements  of  the  constituents  and  not  as 
the  average  of  the  microscopic  stresses.  Mokadam  [100]  following  Guggenheim  [68] 
regarded  the  fluid  pressure  to  be  the  thermodynamic  property  such  that 
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7r  d  V  =  dW  (2.97) 

where  dV  is  the  differential  change  in  the  intrinsic  fluid  volume  and  dW  the  reversi¬ 
ble  work  of  the  fluid  phase. 


Garg  (491  Morland  [103-105]  Pecker  [119]  and  Carroll  [27]  among  others,  intro¬ 


duced  the  notion  of  intrinsic  stresses  for  each  constituent  leading  to 


Ar 


(2.98) 


Tsien  [166]  divided  the  total  stress  into  stress  deviation  and  a  hydrostatic  compo¬ 
nent.  The  hydrostatic  stress  was  expected  to  be  distributed  over  the  solid  and  the  fluid 
in  proportion  to  their  volume  fractions  and  the  solid  was  expected  to  take  the  entire 
stress  deviation.  Terzaghi  [161,1621  as  well  as  Green  [61-67]  assumed  the  partial  stress¬ 
es  to  act  over  the  entire  area  of  any  surface  element.  Further,  assuming  partial  stresses 
to  be  additive,  the  total  stress  in  a  saturated  soil,  assuming  isotropic  fluid  stress,  is 

•u-4'H8  <*»> 

Biot  regarded  the  partial  soil  stress  to  be  the  bulk  stress  acting  over  the  entire 
area  of  internal  surfaces.  In  his  earlier  work  [ill  there  was  no  reference  to  the  area 
over  which  the  fluid  pressure  acts.  Later  [121  Biot  assumed  the  fluid  pressure  to  act 
only  over  the  pore  area.  This  corresponds  to  the  notion  of  the  fluid  pressure  being  an 
intrinsic  quantity.  Thus 

t,  =  ,l''+n(Xr  (2.100) 

■j  'J  'j 


a).  Diffusive  Resistance. 

Diffusive  force,  identified  as  the  interaction  force  by  Truesdell  [163,164]  was 
called  "diffusive  resistance"  by  Green  [62-67]  and  Crochet  [37].  Biot  [18]  described  it 
as  "disequilibrium  force".  The  words  "friction”  and  "drag”  have  also  been  used.  Max- 
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well  [97]  defined  it  as  a  pair  of  equal  and  opposite  forces  acting  on  the  two  constitu¬ 
ents  in  a  binary  mixture.  For  non -chemically  reacting  continue,  in  the  absence  of  iner¬ 
tia  effects,  the  equilibrium  equations  for  the  fluid-saturated  solid  can  be  re-written  as 

$+*(T  -  -t(^-P(2X2)  (2.101) 

Each  side  of  the  equality  represents  interaction  between  the  constituents  and  is  set 
equal  to  components  of  the  diffusive  resistance  vector.  A  set  of  single-constituent  stress¬ 
es  in  equilibrium  can  be  added  to  the  stresses  on  either  side  without  affecting  the  def¬ 
inition  of  diffusive  resistance.  For  hydrostatic  fluid  stress,  the  diffusive  resistance  is 


D,  =  A!  vs 


i<i> 

J.j 


(2.102) 


where 


(t^’+t^  =  0  (2.103) 

the  superposed  bar  indicating  a  set  of  stresses  in  equilibrium.  Maxwell  assumed  the 
interaction  force  to  be  proportional  to  the  densities  of  the  constituents  and  to  the  rela¬ 
tive  velocity.  This  is  easily  seen  to  be  a  special  case  of  the  more  general  relationship 
indicated  by  the  above  definition  of  diffusive  resistance. 


If  inertia  effects  are  included,  for  a  binary  mixture,  p,co  the  quantity  conjugate  to 
the  relative  velocity  in  the  energy  equality  is 

-  p,  =  p !‘X  ff-  bl")  +  pc<!>v“-  (2.104) 

This  definition  is  somewhat  more  general  than  the  one  used  by  Green  [64,65]  where 
the  term  involving  mass  supply  was  ignored.  In  [66],  however,  following  Mills  [99],  a 
general  form  was  stated  for  a  mixture  of  n-1  fluids  and  a  solid. 


Truesdell  [164]  proposed  a  mechanical  theory  of  diffusion  and  showed  that  the 
kinetic  theory,  the  hydrodynamical  theory  and  the  thermodynamical  theory  were  all 
specializations  of  his  general  theory.  According  to  Truesdell;  "diffusion,  being  a  change 
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of  motion,  arises  from  forces;  the  motions  produced  by  these  forces  must  conform  to 
the  principle  of  linear  momentum  applied  to  each  constituent  and  to  the  whole  mix¬ 
ture”.  The  sum  of  supplies  of  partial  momenta  to  the  constituents  must  vanish.  The 
simplest  constitutive  equation  for  supply  of  momentum  would  be  a  linear  dependence 
of  the  partial  momentum  upon  the  quantity  "corresponding”  to  it  in  the  energy  bal¬ 
ance  equation,  viz,  the  relative  velocity  of  the  constituent  with  respect  to  the  others. 
The  restriction  that  the  sum  of  partial  moment  supplies  must  vanish  places  a  restric¬ 
tion  upon  the  coefficients.  Truesdell  showed  that  this  restriction  leads  to  the  necessary 
and  sufficient  condition 

£(1^-1$  =  0  (2.105) 

ri 

where  the  scalar  coefficients  relating  velocity  of  the  kth  constituent  with  its  par¬ 
tial  momentum  are  uniquely  defined  for  k  ^  j  and  must  vanish  for  k  =  j.  For  a  binary 
mixture,  this  implies  symmetry  with  respect  to  j  and  k. 


Mokadam  [100]  proposed  that  constitutive  equations  for  the  diffusive  force  have 
the  form 

p  =  —  JL  v. - T  .  (2.106) 

Crochet  [37l  assuming  the  existence  of  an  energy  function  for  the  mixture, 
showed  that  under  isothermal  conditions  and  in  the  absence  of  chemical  reactions,  the 
constitutive  relations  will  involve  deformation  gradients,  deformation  rates,  velocities, 
and  relative  vorticities.  For  the  linear  theory  of  irrotational  relative  motion,  in  the 
absence  of  chemical  reaction  and  inertia  effects,  and  non-Newtonian  behavior  in  a  bina¬ 
ry  mixture  of  a  fluid  and  a  solid  this  immediately  leads  to  an  expression  of  the  type 

,  (2).  (2)  r, 

D  =  7T  +  p  o  =:  — C  V 

.i  r  >  Pi 


(2.107) 
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The  inverse  form  of  this  equation  is  the  well-known  d'Arcy  flow  rule  [120,141].  Biot 
[11,14-21]  and  earlier  investigators  had  used  this  rule  as  the  starting  point  for  their 
theories.  The  observation  of  linear  dependence  of  one-dimensional  water  flow  and  the 
potential  gradient  has  been  extended  to  two  and  three-dimensional  flow.  The  constant 
of  proportionality  was  enlarged  to  have  the  nature  of  a  second  rank  tensor  transform¬ 
ing  the  potential  gradient  to  the  flux  vector.  The  permeability  tensor  is  generally 
assumed  to  be  symmetric  cone  ponding  to  Biot's  assumption  of  the  existence  of  a  dissi¬ 
pation  function.  Nonhomogeneity  of  the  solid  and  spatial  variations  in  fluid  properties 
have  been  allowed  for  by  aiming  the  components  of  the  permeability  tensor  to  be 
spatially  varying.  The  tensorial  character  has  been  used  to  admit  hydraulic  anisotropy. 
In  case  where  the  solid  as  well  as  the  fluid  are  in  motion,  d'Arcy's  law  has  been 
applied  to  the  relative  velocity  of  the  fluid  with  respect  to  the  solid  matrix  (e.g. 
[14,55]). 
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where  h,,  m,,  v,,  D,,  denote,  respectively,  components  of  the  heat  flux,  the  diffusive 

flux,the  mass  velocity,  and  the  diffusive  force  vector.  fim  is  the  chemical  potential  of 

the  mth  constituent  For  isothermal  flow,  in  the  absence  of  chemical  reactions,  this 
equation  merely  indicates  the  temperature  dependence  of  the  permeability  tensor. 

The  d'Arcy  fluid  flow  equation  has  been  generalized  still  further  [58,591  and 
indeed  forms  part  of  the  general  phenomenological  equations  given  by  Onsager  [118] 

These  express  the  effect  of  simultaneous  presence  of  fields  of  mechanical  pressure,  elec¬ 

tric  potential,  temperature  and  chemical  concentration.  The  general  relationship  is 
expressed  as 

J,  =  L/j  (2.110) 

where  J,  are  the  fluxes  viz^  mass  flow,  heat  flow,  electric  current,  chemical  diffusion. 
Fj  are  the  Prigogine  forces  [59,1181  L,  are  components  of  a  positive  definite  symme¬ 
tric  tensor.  The  components  L#  have  to  satisfy  the  Curie-Prigogine  principle  and  in 
some  cases  of  system  symmetry,  L,  would  vanish  where  the  tensorial  rank  of  the 
"forces”  Fj  and  the  "fluxes”  J,  are  not  the  same.  Onsager  [118]  expected  the  relationship 
to  be  symmetric  Le,  L,j  =  Lr  Evidently,  the  Onsager  equation,  proposed  originally  for 
small  perturbations  on  an  equilibrium  state,  is  a  restricted  type  of  relationship,  assum¬ 
ing  a  quadratic  form  for  the  entropy  rate  function.  In  general,  J,  can  be  treated  as 
functionals,  and  may  depend  linearly  or  nonlinearly  on  the  spatial  gradients,  of  all 
orders,  of  the  potential  fields  and  their  history.  This  has  been  discussed  by  Coleman 
[36]  in  the  case  of  heat  conduction.  The  same  type  of  reasoning  would  apply  to  other 
flux  phenomena,  and  the  only  thermodynamic  restriction  is  that  the  scalar  product  of 
flux  and  force  be  non-negative  [164]. 

Ghaboussi  [54]  pointed  out  that  Biot's  [17,18]  formulation  for  momentum  balance 
may  be  regarded  as  a  generalization  of  (2.107)  to  include  inertia  effects  in  the  body 
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force  term  giving,  for  no  chemical  reaction  and  isotropic  fluid  stress,  a  specialization  of 
the  more  general  relationship, 

p,  =  TTj  +  pt2)(  bj2>— fj2>)  =  —  CjjVj  (2.111) 

Using  the  definition  of  p,,  (2.69),  a  general  relationship  based  on  the  correspondence  of 


p}"  and  the  relative  velocity  v[°— vf5  is 


(2.112) 


b).  Stresses 

It  has  been  difficult  to  define  the  dependence  of  components  of  the  partial  stress 
tensors  upon  the  kinematic  variables  and  densities  of  the  constituents  of  a  mixture.  In 
his  earlier  work,  Adkins  [1,2]  assumed  that  the  stress  in  each  constituent  depended 

upon  the  density  and  the  kinematic  quantities  associated  with  that  constituent  only. 
Green  [60]  admitted  interdependence  of  partial  stresses  of  each  constituent  upon  the 
kinematics  of  all  This  was  in  line  with  the  principle  of  equi presence.  However,  to 

make  the  independent  variables  distinct  for  various  constituents.  Green  [60]  stipulated 
that  the  partial  stresses  for  any  constituent  would  depend  upon  the  densities,  velocities 
and  antisymmetric  deformation  of  all  constituents  but  only  on  the  symmetric  part  of 
the  deformation  gradient  of  the  constituent  itself.  Invariance  of  stress  under  superposed 
rigid  body  motions  and  under  superposed  uniform  rigid  body  angular  velocities  of  the 
mixture  as  a  whole  showed  that,  for  a  binary  mixture,  the  velocities  and  the  rotation 
tensors  must  occur  as  difference  terms  in  the  set  of  independent  variables.  Mokadam 
[100]  assumed  the  stress  tensor  to  be  a  linear  function  of  velocity  gradients.  Carroll 

[27]  and  Morland  [103-106],  among  others,  also  assumed  that  the  shear  traction  is  car¬ 
ried  by  the  bulk  solid  and  is  related  to  changes  in  pore  geometry.  The  hydrostatic 

stress  components  in  the  solid  and  the  fluid  were  expected  to  cause  intrinsic  volume 
changes  in  the  constituents. 
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Bedford  [9]  pointed  out  that  fluid-saturated  porous  media  fall  in  the  class  of 
immiscible  mixtures.  The  constituents  of  such  mixtures  remain  physically  separate  on 
a  scale  which  is  large  in  comparison  with  molecular  dimensions.  This  immiscibility 
has  two  important  consequences.  Because  of  physical  separation,  in  some  local  sense, 
each  constituent  will  obey  the  constitutive  relations  for  that  constituent  alone.  Also, 
the  constituents  intrinsically  have  microstructure  defined  by  the  inter-faces  which  sepa¬ 
rate  the  constituents.  To  set  up  macroscopic  constitutive  relations,  one  approach  would 
be  to  postulate  these  relations  directly  as  described  above.  Another  alternative  would  be 
to  relate  macroscopic  behavior  to  intrinsic  properties  of  the  constituents.  The  simplest 
theories  involve  the  volume  fraction  of  the  constituents  in  addition  to  the  usual  vari¬ 
ables.  Morland  [103]  pointed  out  the  meaning  of  deformation  and  stress  associated  with 
the  continuum  model  of  each  constituent.  In  particular,  the  partial  density  variation  is 
not  the  density  variation  of  the  constituent  since  the  mixture  postulate  eliminates  refer¬ 
ence  to  the  actual  volume  occupied  by  each  constituent  in  an  immiscible  mixture. 


Terzaghi  [161]  introduced  the  concept  of  effective  stress.  This  was  defined  to  be 
the  stress  component  causing  deformations  of  the  soil.  For  hydrostatic  pore-water  pres¬ 
sure  acting  on  incompressible  soil  grains,  the  entire  deformation  of  the  soil  was 
assumed  to  be  due  to  changes  in  the  pore  volume  and  pore  geometry.  For  this  case, 
TeTzaghi  [161]  called  the  partial  solid  stress  the  effective  stress  related  to  deformation 
of  the  solid.  Biot  [ll]  regarded  the  total  stress  and  the  fluid  pressure  as  the  mechani¬ 
cal  variables.  It  was  found  (e.g.  [l  1 3D  that  the  fluid  pressure  did  in  fact  influence  the 
effective  stress-strain  relationship  when  the  solid  grains  had  compressibility  comparable 
to  that  of  the  matrix  as  a  whole  and  the  fluid  was  not  incompressible.  To  allow  for 
this  the  effective  stress  related  to  deformation  was  defined  as 

(2.113) 

=  t(1)  +  (  1  —  chrS 


t  =  t.  —  cwS 

>j  >j  u 


(2.114) 
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where  c  -  1  implies  Terzaghi's  definition  and  c  -  0  would  correspond  to  total  stress 
being  regarded  as  effective.  Suklje  [160]  discussed  selection  of  appropriate  values  of  c. 
Schiffman  [143]  expected  c  to  be  between  the  value  of  porosity  and  1.  This  was  based 
on  the  assumption  that  the  pore  fluid  pressure  may  not  act  over  the  entire  area  of 
surface  elements  but  only  over  a  part.  This  fractional  area  of  action  of  the  fluid  pres¬ 
sure  would  be  bounded  below  by  the  porosity  and  above  by  1.  Several  investigators 
confuse  the  effective  stress  with  the  partial  stress.  This  is  due  to  the  dual  definition 
originally  given  by  Terzaghi  The  term  "effective  stress”  used  in  this  report  is  the 
stress  component  causing  deformations  of  the  solid  and  is  thus  defined  completely  by 
these  deformations.  The  partial  stress  in  the  solid  could  conceivably  be  related  to  quan¬ 
tities  other  than  the  deformations  of  the  solid.  Verruijt  [167,168]  used  the  term 
"intergranular  stress”  for  the  difference  between  the  total  stress  and  the  intrinsic  pore- 
water  pressure  assumed  to  act  over  100  percent  area.  The  relationship  of  "intergranular 
stress”  and  the  "effective  sress”  with  the  partial  stress  which  appears  in  equations  of 
balance  of  momentum  needs  to  be  established.  Considering  the  solid  grains  to  be  com¬ 
pressible,  Nur  [113]  derived  the  equation 

c  =  1  -  K/Ki  (2.115) 

where  K,  K,  are,  respectively,  the  bulk  and  the  intrinsic  compressibility  of  the  solid. 
For  incompressible  grains  and  highly  deformable  pore  space,  K/K,=0,  and  Terzaghi's  two 
definitions  coincide.  Suklje  [160]  suggested,  without  proof, 

c  =  l-n(1>-£-  (2.116) 

lx 

Schiffman  [143]  gave  a  more  general  form  allowing  fluid  pressure  to  be  a  second  rank 
tensor  and  c  a  fourth  rank  tensor,  i.e. 


t.  =  t.  -  A....t(*) 

ij  ij  ilij  kl 


(2.117) 


Carroll  [27]  carried  out  a  similar  development.  These  approaches  were  based  on  superpo¬ 
sition  of  effects  of  the  hydrostatic  stress  and  the  stress  deviation.  Carroll  [27]  deter¬ 
mined,  for  the  linear  case,  the  relation 

where  we  components  of  the  elasticity  tensor  for  the  dry  solid  material  and 
those  of  the  intrinsic  compliance  under  hydrostatic  stress.  Biot  [18]  defined  effective 
stress  as 

t.  =  V-oirV  (2.119) 

ij  ij  u 

and  assumed  t'(J  to  be  the  quantity  related  to  solid  deformation.  This  coincides  with 
Terzaghi’s  concept  of  effective  stress  for  a  =  1. 

Garg  [491  following  Haimson  [73l  proposed  a  dual  definition  for  effective  stress. 
For  strength  of  rock,  he  would  set  c  *  1  but  for  constitutive  relations  another  value 
of  c  would  be  used. 

Carroll  [28]  introduced  intrinsic  solid  stress  on  the  solid  particles  and  an  effective 
stress  influencing  deformation  of  the  pore  space.  Kenyon  [83,84]  also  considered  the 
effect  of  grain  and  fluid  compressibilities  and  introduced  material  parameters  to  charac¬ 
terize  this  dependence.  Contact  stress  in  the  solid  and  the  bulk  stress  independent  of  K, 
were  used. 

For  large  deformation,  an  incremental  form  of  the  stress  tensor  was  introduced 
by  Biot  [20].  Carter  [30]  and  Prevost  [121]  used  the  Jaumann  stress-rate  to  ensure 


frame  indifference. 
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Form  of  the  Stress-Strain  Relations. 


Gibson  [55]  expected  the  effective  stress  to  depend  upon  the  deformation  or  the 
rate  and  history  of  deformation  of  the  solid  skeleton.  The  fluid  pressure  was  expected 
to  depend  upon  the  fluid  density  in  an  isothermal  system. 

Tsien  [167]  proposed  a  linear  elastic  isotropic  relation  for  the  partial  soil  stress  in 
terms  of  the  soil  strain  using  Terzaghi's  definition  Le.  c  =  l  in  (2.113).  Biot  [l  l] 
assumed  a  quadratic  energy  function  in  the  change  in  water  content  per  unit  volume 
of  the  solid,  and  the  soil  strains  leading,  for  isotropic  linear  elastic  soil  and  incompres¬ 


sible  fluid,  to 


n  =  Metk  +  N0 


*U  =  2Meij  +  Xekk8ij  +  Me61. 


(2.120) 


(2.121) 


In  later  work,  [17],  the  total  stress  was  replaced  by  the  partial  solid  stress.  In  exten¬ 
sion  to  anisotropic  elastic  [15,16]  solid  and  compressible  fluid  the  relationships  for  the 


total  stress  and  the  intrinsic  pore-water  prssure  were  stated  as 


(2.122) 


(2.123) 


An  alternative  form  of  the  above  relationships,  assuming  My  =  aMS^  and  ty  =  t'(J  +  crrrS^ 


[16,18,191  is 


tij  =  Euifu  +  oMSjj^u+C) 


it  =  aMe(kk  +  M< 


(2.124) 


(2.125) 


Here  a  is  a  measure  of  compressibility  of  the  solid  particles.  This  form  was  used  by 
Ghaboussi  [53,54]  for  development  of  finite  element  solution  procedures.  Garg  [50]  wrote 
the  constitutive  relations  for  an  isotropic  system  in  the  form 
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C  =  aeu5y  +  ce2S.i  + 


(2.126) 


_  -  „<*V  -  _  JD  ,  .(2) 

ir  -  n  ir  -  ce^+be^ 


Here  y.  is  the  bulk  shear  modulus  of  the  porous  solid  and  a,b,c  are  functions  of  the 
volume  fractions  of  the  constituents,  the  bulk  modulus  of  the  porous  solid,  and  the 
intrinsic  bulk  moduli  of  the  fluid  and  the  solid.  The  constants,  for  isotropic  elasticity 
have  been  shown  to  correspond  to  Biot's  and  to  depend  upon  the  properties  of  the  con¬ 
stituents  and  their  volume  fractions. 


Similar  construction  was  used  for  viscoelastic  soils  [12].  In  [20]  the  fluid  strain 
was  again  replaced  by  the  change  in  water  content.  The  same  concept  was  extended  to 
the  case  of  finite  elastic  deformation  [21]. 


Lubinski  [93]  assumed  that  the  total  strain  of  a  bulk  porous  solid  can  be 
expressed  as  a  summation  of  the  strains  due  to  pore-water  pressure  and  strains  due  to 
stress  acting  on  the  solid  skeleton.  He  proposed  relations  of  the  type 

JO  —  c  J 0  f  J *) _ ..  (2127) 


All  _  C  ,  /  Hi  -S 

Vi  =  Etiyeti  +(n  -V)ir8ij 


n  =  Me^1 1  +  N  ( n(2)-  n(Q2))  (2.1 28) 

where  y,  M,  N,  are  material  constants.  Krause  [86j  added  terms  to  the  right  side  of 
(2.120)  to  reflect  linear  dependence  of  the  fluid  pressure  on  the  deformation  rate  of 
the  fluid.  This  assumes  a  viscous  component  for  fluid  flow.  Adkins,  in  his  earlier 
theory  [1,21  assumed  that  the  stress  in  each  constituent  depended  only  on  the  density 
and  the  kinematic  quantities  associated  with  only  that  constituent.  Nur  [113]  assumed 
the  effective  stress,  given  by  (2.113)  along  with  (2.115),  to  be  related  to  the  soil 
strains.  This  admitted  a  certain  dependence  of  the  partial  soil  stress  upon  the  fluid 
pressure.  Explicitly, 
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(l) 


lij  =  tlj  1  C^ij  = 


(2.129) 


Hence 


^C^-d-c)^  (2.130) 

These  approaches  were  based  on  the  superposition  of  effects  of  the  hydrostatic  stress 
and  the  shear  stress.  Carroll  [27]  determined,  that  for  the  linear  case,  (2.130)  would  be 

<!’  -  ■wff-'wSl'  (2-131) 

where  are  components  of  the  elasticity  tensor  for  the  dry  9olid  material  and 
are  those  of  the  intrinsic  compliance  under  hydrostatic  stress.  This  formulation  was  spe¬ 
cialized  to  allow  for  the  presence  of  internal  symmetries.  For  isotropy,  the  formulation 
reduces  to  Nut's  [113L  Schiffman  [143]  proposed  a  more  general  form  of  (2.131)  viz, 

■!,“  -  Burf’ -(»„*, -^4*  (2-132) 

The  quantity  A,^,  was  termed  the  soil-water  interaction  tensor.  Garg  [49]  obtained  a 
relationship  between  the  intrinsic  and  the  bulk  behavior  of  rocks  under  hydrostatic 
stress. 


In  extending  the  theory  to  the  nonlinear  case,  Westmann  [170]  assumed  the  par¬ 
tial  solid  stress  to  be  a  function  of  the  deformation  tensor  for  the  solid  and  the  rate 
of  deformation  (Eulerian  description)  of  the  fluid.  The  fluid  stress  was  expected  to  con¬ 
sist  of  a  hydrostatic  component  and  another  component  depending  upon  the  same  quan¬ 
tities  as  the  partial  solid  stress.  It  was  noted  that  in  this  formulation  it  would  be  dif¬ 
ficult  to  design  experiments  to  evaluate  the  parameters.  A  simplification  proposed 
assumed  the  fluid  pressure  to  be  hydrostatic  and  related  to  the  velocity  field  through 
d'Arcy's  law.  This  is  similar  to  Sandhu's  [131-133]  argument  that  the  constitutive  equa¬ 
tion  for  diffusive  resistance  is  a  sufficient  relationship  between  the  fluid  partial  stress 
and  kinematics  of  the  mixture.  Westmann  [170]  wrote  relative  velocity  as  a  function 
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of  the  partial  fluid  stress,  the  diffusive  resistance  and  the  Cauchy  deformation  tensor 
for  the  solid.  This  would  reflect,  among  other  factors,  the  dependence  of  permeability 
on  the  porosity  of  the  solid. 


Morland  [103]  and  Garg  [49]  assumed  the  intrinsic  stress  in  each  component  to  be 
a  function  of  the  deformation  of  that  constituent  only  and  having  the  same  form  as 
for  a  single  material.  Thus 

t(k)*  =  f  [e(k)*]  (2.133) 

jj  mm 

The  bulk  stresses  and  deformations  were  expected  to  be  related  to  the  corresponding 
intrinsic  quantities  through  scaling  functions.  Thus,  the  bulk  stress  in  the  kth  constit¬ 
uent  is 


(k)  _  (k)Jk)* 

t..  —  n  i.. 
jj  jj 


(2.134) 


For  linear  isotropic  elastic  rock 
JD  _  (iV,( HMD* 


,(D  _  dUD 

x  ii  ~  n 


'ii 


(2.135) 


where  the  subscript  D  denotes  "dry”  rock-  The  bulk  deformation  gradient  was  related 
to  the  effective  deformation  gradient  as: 

f‘.  »  [  a(IW0l),/3  (2.136) 

The  relation  between  the  deformation  gradient  and  the  partial  soil  stress  was  expected 
to  be  the  same  function  f  as  for  the  intrinsic  quantities.  The  isotropic  pressures  in  the 
fluid  and  the  solid  were  assumed  to  depend  upon  the  volumetric  strain  of  both  the 
constituents,  i-e. 

t..  *  ae<l)  +  be(2)  (2.137) 

jj  mm  mm 


(1)  ,  .  (2) 
ir  =  ce  +  de 

mm  mm 


(2.138) 


However,  unlike  Biot,  the  existence  of  an  energy  function  was  not  postulated  so  that 
the  constants  b  and  c  in  (2.137)  and  (2.138)  do  not  have  to  be  equal.  Morland 
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l 

I 
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expressed  the  coefficients  in  terms  of  compressibilities  of  the  fluid  and  the  solid,  and 
the  bulk  shear  modulus  of  the  mixture. 

« 

Carroll  [29]  postulated  the  following  relations  for  a  fluid-saturated  solid: 


L  Relation  between  stress  in  the  mixture  and  stresses  in  the  constituents: 

■jta  =  in^tJWV  (2.139) 


ii.  Solid  stress-strain  law 

1  it)  _  _v  aV0  • 

where  the  symbol  A  indicates  change  in  the  quantity  following  it. 

iii.  Effective  stress-strain  law: 


« 

(2.140) 


« 


it  itu  = 


(2) 


3  ii  3  ii 


(l) 


(2.141) 


Combining  these  relationships  they  obtained  the  bulk  relations  for  the  mixture. 


Thermodynamic  considerations. 


i 


Adkins  [3]  and  Green  [60]  admitted  interdependence  of  stress  of  each  constituent 
upon  the  kinematics  of  alL  This  was  in  line  with  the  principle  of  equi  presence  stated 
by  Truesdell  [164].  In  application  to  elastic  materials,  the  existence  of  an  energy  func¬ 
tion  for  the  mixture  was  assumed  by  Biot  [11,14-16,19-21].  This  has  been  consistently 
followed  by  numerous  investigators  (e.g.  [6, 7,25, 37,39, 62-67, 159D-  Sandhu  [131-134] 
pointed  out  that  as  the  mixture  could  not  be  regarded  as  a  continuum  in  motion,  it 


< 


was  inappropriate  to  assume  energy  functions  for  it  in  the  form  that  has  been  popular. 


Sandhu  [131-134],  Westmann  [171]  and  Morland  [103-105]  followed  Adkins'  [2] 
original  idea  that  the  stresses  in  each  constituent  depend  upon  the  kinematics  of  only 
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that  constituent.  However,  Moriand  [103-105]  used  this  for  the  intrinsic  rather  than 
the  bulk  stresses.  This  brings  back  some  dependence  of  the  partial  solid  stress  upon  the 
fluid  pressure  because  the  porosity  was  postulated  to  be  a  linear  function  of  the  par¬ 
tial  stresses. 

Seeking  constitutive  equations  for  internal  energy,  entropy,  heat  flux  vector,  par¬ 
tial  stresses  and  diffusive  resistance  vectors  p,  and  the  quantity  q— h^  in  the  case  of 
mixture  of  two  Newtonian  fluids,  Green  [62]  assumed  these  to  depend  upon  the  densi¬ 
ties  of  the  constituents,  the  velocities,  the  gradients  of  velocities,  and  the  temperature. 
For  the  heat  flux  vector,  the  temperature  gradient  replaced  the  deformation  gradients. 

Discussing  thermodynamics  of  fluid  flow  in  a  rigid  porous  medium,  Mokadam 
[101,102]  pointed  out  that  d'Arcys  law  is  valid  only  for  isothermal  flow  in  which  the 
inertial  and  viscous  effects  are  negligible.  Also  that  the  rate  of  entropy  production 
must  be  non-negative  separately  for  terms  involving  quantities  of  different  tensorial 
ranks.  Crochet  [37]  applied  thermodynmic  considerations  to  the  flow  of  a  fluid  through 
an  elastic  solid.  Atkin  [6]  explicitly  stated  the  form  of  these  constitutive  assumptions 
for  flow  of  a  fluid  through  an  elastic  solid.  He  also  presented  an  alternative  method 
of  deriving  the  linearized  theory  of  elastic  solid-viscous  fluid  mixtures  and  the  thermo¬ 
dynamic  restrictions  imposed  on  this  theory  by  the  entropy  production  inequality.  In 
later  work.  Green  [67]  based  the  thermodynamic  restrictions  on  the  behavior  of  each 
constituent  on  the  requirement  that  suitable  combinations  of  the  equations  for  individu¬ 
al  constituents  should  yield  a  single  entropy  production  inequality  for  the  mixture  as  a 
whole.  Bowen  [23]  noted  that  these  formulations  lead  to  the  result  that,  in  equilibri¬ 
um,  the  partial  free  energy  density  of  a  given  constituent  is  independent  of  the  defor¬ 
mations  of  the  other  constituents.  Also  that  such  independence  fails  to  be  confirmed  by 
experiments  on  fluid  mixtures.  Muller  [108]  showed  that  if  gradients  of  densities  of 
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the  constituents  were  included  among  the  constitutive  variables,  this  difficulty  would 
not  arise.  Some  investigators  have  proposed  use  of  an  entropy  production  inequality  for 
each  constituent.  Bowen  [23]  considers  this  to  be  too  special 

Morland  [103-105]  did  not  assume  the  existence  of  an  energy  function  for  the 
mixture  but  still  admitted  interdependence.  This  implies  admitting  a  possibly  nonsy in- 
met  ric  constitutive  relationship  of  the  type  proposed  earlier  by  Schiffman  [143].  For  c 
-  porosity,  the  constitutive  equations  for  stresses  become  uncoupled. 

Crochet  [37]  started  by  admitting  fairly  general  constituive  assumptions  in  line 
with  Truesdell's  [163]  principle  of  equipresence,  and  then  determined  the  restrictions 
placed  upon  these  general  constitutive  relations  by  thermodynamic  considerations.  This 
approach  is  similar  to  that  used  by  Noll  [ill]  and  Coleman  [36].  Crochet  [37]  showed 
that  the  restriction  of  nonnegative  entropy  production  requires  that  the  entropy  and  the 
internal  energy  be  independent  of  the  deformation  rates,  relative  velocity  and  the  temp¬ 
erature  gradient  Green  [66]  extended  Crochet's  [37]  work  to  anisotropic  solids  and  to 
include  initial  stresses. 

d.  Constitutive  Relations  for  Porosity. 

Gibson  [55]  treated  porosity  as  a  function  of  effective  stress  and  proposed  compliance 
relationships  in  the  form 

n(2)  =  fCt(J,X.,t)  (2.142) 

Walsh  [169]  regarded  the  pore  space  or  the  volume  fraction  of  the  pores  to  be  a 
function  of  the  solid  stress.  This  leads  to  the  relationship; 

dn,2>  1  dV*2)  n(2)  dV  1  dV2)  V2)  dV 

dt(,)  V  dt<»  V  dt<»  "  V  dtC)  v2  dt(1) 

U  il  U  il  JJ 


(2.143) 
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Garg  [49]  pointed  out  that  Walsh's  analysis  was  acceptable  for  very  dense  rock.  For  the 
general  case,  they  introduced  bulk  and  intrinsic  solid  compressibilities  and  set  up  more 
general  expressions.  They  proposed  constitutive  equations  for  porosity  based  on  the  exis¬ 
tence  of  an  energy  function. 

Aif antis  [4]  assumed  effects  of  changes  in  fluid  pressure  and  the  solid  stress  to 
be  additive  and  proposed  a  compliance  relationship 

An(2)  =  aAir  +  bAt^  (2.144) 

Assuming  the  intrinsic  properties  of  any  constituent  are  not  affected  by  the  pres¬ 
ence  of  the  other  constituents,  Morland  [103]  proposed  a  constitutive  equation  for  poros¬ 
ity  in  the  form 

n(2)  =  n(02)(  1  +  at'!1*  bir  )  (2.145) 

To  include  dilatancy,  the  relationship  was  generalized  further  to 

n<2)  =  n(02)(  1  +  a.-tjj’  +  brr  )  (2.146) 

2-3 .2.7  Comments. 

Various  approaches  to  description  of  the  constituents  and  the  mixture  as  well  as 
description  of  their  motion,  formulation  of  the  equations  of  balance  of  mass,  linear 
momentum,  angular  momentum,  and  rate  of  energy,  and  the  constitutive  relations  have 
been  discussed.  In  most  theories  of  mixtures,  deformation  is  referred  to  an  initial  con¬ 
figuration  for  each  constituent  and  motion  to  the  place  coordinates.  Also  the  equations 
of  balance  are  written  for  a  fixed  volume  in  space.  This  approach  may  not  be  conven¬ 
ient  for  soil- water  mixtures. 

Truesdell  [164,166]  postulated  equations  of  balance  of  mass,  linear  momentum, 
moment  of  momentum  and  energy  such  that  the  form  of  the  equations  was  the  same 
for  each  constituent  and  for  the  mixture.  The  notion  of  motion  of  a  mixture  as  a 
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whole  was  introduced.  Indeed,  Truesdeli  would  require  the  form  of  the  equations  of 
balance  for  the  mixture  to  be  the  same  as  for  a  single  material  To  accomplish  this 
identity  of  form,  the  total  stress  tensor,  the  total  heat  flux  vector,  and  the  specific 
energy  supply  had  to  be  specially  defined  and  did  not  equal  the  stun  of  the  corre¬ 
sponding  quantities  for  the  constituents.  The  specific  energy  (internal  plus  kinetic)  of 
the  mixture  was,  however,  equal  to  the  sum  of  the  corresponding  quantities  associated 
with  the  constituents.  The  existence  of  the  mixture  as  a  continuum  in  motion  with 
acceleration  derived  from  the  barycentric  velocity  is  implied  in  this  line  of  thought. 
The  analysis  was  founded  on  the  so-called  fundamental  identity  involving  "material 
derivatives  of  the  mean  value”.  Whereas  these  can  be  accepted  as  hypothetical  entities 
for  simplification  of  analysis,  it  is  difficult  to  assign  a  physical  meaning  to  them.  This 
material  rate  has,  in  Atkin's  [7]  words,  "no  particular  physical  significance".  This  is  so 
because  the  rate  is  executed  not  on  a  material  particle  but  on  a  center  of  mass.  The 
mixture,  at  any  instant  of  time,  has  been  constructed  by  superposition  of  constituent 
particles,  is  not  a  set  of  particles,  consists  only  of  centers  of  mass,  and  is  defined  only 
for  the  particular  instant  of  time.  It  cannot  be  regarded  as  a  continuum,  consisting  of 
a  set  of  non-penetrating  particles,  in  motion.  Atkin  [7]  pointed  out  that  the  mixture 
density  cannot  be  associated  with  a  material  in  the  physical  sense.  Sandhu  [133,134] 
pointed  out  that  the  mixture  defined  above  has  a  physical  existence  only  in  the  case 
of  no  relative  motion  between  constituents.  For  this  special  case  the  mixture  will  have 
motion  and  deformation  as  a  material  body  and  the  development  of  equations  of 
motion  for  the  mixture  is  meaningful,  for  example,  in  the  post-liquefaction  stage.  How¬ 
ever,  in  study  of  wave  propagation  leading  to  liquefaction  of  soils,  it  is  of  little  inter¬ 
est.  If  relative  motion  is  present,  the  mixture  does  not  satisfy  the  axiom  of  continuity 
and  its  corollary,  the  principle  of  impenetrability.  Accordingly,  the  mixture  density, 
momentum,  moment  of  momentum  and  energy  defined  by  Truesdeli  are  only  mathe- 
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matical  entities  without  any  physical  interpretation.  These  quantities  cannot  be  regarded 
as  functions  of  time  associated  with  a  set  of  physical  particles.  Truesdell's  third  postu¬ 
late,  therefore,  appears  to  be  irrelevant  to  the  development  of  a  general  theory. 

Chao  [34]  noted  that  combining  the  balance  equations  of  constituents  to  obtain  the 
balance  equations  of  the  mixture  can  lead  to  errors.  As  an  example,  the  absence  of 
inertial  coupling  forces  in  the  momentum  equations  of  the  constituents  was  cited.  This 
is  contrary  to  Biot's  "mass  coupling"  assumption. 

On  the  other  hand.  Green  [62]  considered  the  concepts  of  stress,  heat  flux,  and 
energy  supply  to  be  primitive  to  each  constituent  and  to  the  mixture  as  a  whole  as 
well.  The  additive  property  of  stress,  heat  flux,  and  the  energy  supply  was  postulated 
and  the  balance  laws  derived  from  the  frame  invariance  of  a  rate  of  energy  equality. 

The  energy  density  of  the  mixture  was  seen  to  be  different  from  the  sum  of  energy 

densities  of  the  constituents.  This  was  attributed  to  interaction  between  the  constituents. 
Green  [66]  established  a  relationship  between  these  quantities. 

The  balance  equations  due  to  Truesdell  [164,166]  and  to  Green  [66]  have  similar 
form  and  are  essentially  equivalent  but  the  quantities  appearing  in  the  two  sets  have 

different  interpretations  based  upon  the  relationships  postulated  between  the  quantities 

associated  with  the  constituents  and  with  the  mixture.  Gurtin  [69,70]  and  Morland 
[103]  support  the  additivity  of  partial  stresses  on  the  ground  that  tractions  are  additive 
and  Cauchy's  stress  principle  should  hold  for  total  stress  and  total  traction  as  well  as 
for  the  constituents.  In  Green's  theory,  the  equations  of  mass  and  momentum  balance 
are  derived  from  the  material  frame  invariance  of  a  rate  of  energy  equality.  In 
another  discussion  the  heat  fluxes  and  the  energy  supply  were  assumed  to  be  additive. 
In  another,  more  recent  version  of  the  theory,  Green  [64]  made  the  role  of  interactions 
between  constituents  explicit  by  writing  the  rate  of  energy  equality  for  each  constitu¬ 


ent. 
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Bowen  [24]  postulated  the  point  form  of  the  rate  of  energy  equality  and  pointed 
out  the  effect  of  certain  approximations.  He  claimed  that  Green's  theory  is  a  special 
case  of  Truesdeil's.  This  is  not  true.  The  interpretation  of  quantities  appearing  in 
Truesdell’s  equations  is  quite  different  from  that  of  similar  quantities  in  Green's  theory 
because  the  two  formulations  are  based  on  different  definitions  for  the  quantities  asso¬ 
ciated  with  the  mixture  in  terms  of  those  for  the  constituents  rather  than  due  to  any 
approximation. 

Westman  [170]  noted  that  while  writing  mass  continuity  relations  care  must  be 
exercised  because  the  volume  of  each  continuum  phase  is  not  the  same  as  the  true  vol¬ 
ume  of  each  material  For  the  case  of  initial  stresses,  equilibrium  must  be  satisfied  in 
the  initial  as  well  as  the  deformed  configuration  realized  after  incrementation  of  stress. 

In  recent  work  by  Gurtin  [69,70]  Oliver  [116,117],  Williams  [171],  and  Sampaio 
[129,130]  the  equations  of  balance  of  momentum  and  energy  differ  from  those  of 
Truesdell  [163,165]  and  Kelly  [82]  They  showed  that  extension  of  the  traditional  theory 
for  single  materials  to  mixtures  by  simply  replacing  the  forces  by  "total  force”  is  inad¬ 
equate  to  express  balance  of  forces  for  other  than  pure  constituents.  In  addition  to  the 
partial  stress  for  each  constituent,  they  obtained  embedding  stresses  governed  by  addi¬ 
tional  balance  of  force  equations. 

Several  investigators  have  introduced  volume  fractions  as  additional  variables  in 
theories  for  compressible  materials.  The  balance  equations  have  been  written  in  terms 
of  relative  motion  and  porosity,  which  is  essentially  a  measure  of  relative  deformation. 
Fukuo's  [46]  equations  of  mass  balance  may  be  regarded  as  an  extension  of  Gibson's  [55] 
approach,  of .  referring  to  the  fixed  set  of  particles  in  the  reference  configuration. 
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In  Gibson's  theory  of  one-dimensional  nonlinear  consolidation,  the  equation  of 
equilibrium  of  vertical  forces  associated  with  reference  volume  of  the  solid  in  the  ref¬ 
erence  configuration  were  written  for  the  current  configuration.  The  density  was  relat¬ 
ed  to  the  of  the  constituents  in  the  current  configuration  which  in  turn 

depended  upon  the  volume  fractions  of  the  constituents.  Beacause,  in  this  representation, 
the  description  of  the  solid  phase  is  unaffected  by  deformation,  the  equation  of  mass 
continuity  for  the  solid  is  simply  the  equation  relating  the  current  density  to  the  ref¬ 
erence  density  of  the  solid. 

The  definition  of  diffusive  resistance  or  the  interaction  force,  given  by  Green 
[62-67]  appears  to  be  appropriate,  as  also  writing  constitutive  relations  for  it  in  terms 
Of  relative  velocities  of  the  constituents.  For  the  simple  case  of  a  binary  mixture,  e.g.  a 
saturated  soil,  the  linear  dependence  of  the  diffusive  resistance  upon  the  relative  veloci¬ 
ty  is  essentially  a  statement  of  the  phenomological  observation  by  d'Arcy.  Biot  assumed 
the  existence  of  a  dissipation  function,  quadratic  in  relative  velocity.  This  corresponds 
to  the  assumption  of  a  linear  dependence  of  velocity  upon  the  pressure  gradients.  Moka- 
rlatn  [100]  pointed  out  that  d'Arcy's  Law  is  valid  only  for  isothermal  flow  in  which 
the  inertial  and  viscous  effects  are  negligible.  Generalizations  of  d'Arcy's  law  to  thermo- 
mechanical  mixtures  and  simultaneous  mechanical  and  chemical  diffusion  along  Onsager's 
principle  form  a  part  of  the  general  theory  of  mechanical  diffusion  presented  by 
Truesdell  [165]. 

For  stresses  in  the  constituents,  several  descriptions  have  been  used.  Constitutive 
equations  need  to  be  written  for  the  partial  stresses  that  appear  in  the  equilibrium 
equations.  Biot  [18,19]  wrote  the  equations  of  momentum  equilibrium  for  the  mixture 
using  the  total  stress.  In  soil  mechanics,  it  is  well  known  that  constitutive  equations 
for  the  total  stress  are  very  sensitive  to  the  pore-water  pressures  and  the  effective 
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stress  is  preferred  for  the  purpose.  Terzaghi  [161]  used  a  dual  definition  for  effective 
stress.  It  appears  reasonable  to  define  it  as  the  component  of  soil  skeleton  stress  which 
is  related  to  the  kinematics  of  the  soil  alone.  According  to  the  other  definition  the 
effective  stress  is  the  difference  between  the  total  stress  and  the  intrinsic  pore-water 
pressure.  Verruijt  [167,168]  calls  it  intergranular  stress.  This  definition  appears  to  be 
unnecessary  except  in  the  case  of  the  so-called  double-porosity  soils  in  which  the  soil 
grain  compressibility  is  taken  care  of  separately  from  the  deformation  of  the  soil  as  a 
whole.  The  relationship  between  the  deformations  of  the  grains  and  the  voids  on  the 
one  hand  and  the  total  soil  mass  on  the  other  has  been  established  for  the  double¬ 
porosity  materials  by  several  investigators  [e.g.  Carroll]. 

Biot's  [18]  and  Crochet's  [37]  assumption  of  the  existence  of  an  energy  function 
for  the  soil-water  mixture  is  open  to  objection.  The  mixture  consists  of  centers  of  mass 
and  not  non-  penetrating  material  particles.  As  such  it  does  not  satisfy  the  axiom  of 
impenetrability  and  it  is  not  correct  to  assign  deformation,  material  rates,  energy  etc.  to 
this  entity.  Thus  the  stress-strain  relationships  for  the  fluid  and  the  solid  may  not  be 
derived  from  an  energy  function.  Biot's  [18,19]  theory  based  upon  the  existence  of  an 
energy  function  quadratic  in  the  strains  of  the  soil  and  the  water  content  or  the  den¬ 
sity  or  the  isotropic  strain  in  the  fluid  is  thus  arbitrary  and  restrictive.  Garg's  [48,49] 
formulation  also  is  similar  and  the  constants  can  be  shown  to  correspond  to  Biot's. 
Garg  related  these  to  the  properties  of  the  constituents  and  the  volume  fractions. 

The  argument  over  whether  the  stresses  in  any  constituent  depend  upon  the  kine¬ 
matics  of  that  constituent  only  or  that  of  all  the  constituents  is  easily  resolved  noting 
that  the  partial  stress  need  not  coincide  with  the  effective  stress.  The  difference 
would  depend  upon  the  compressibility  of  the  fluid  and  the  9oil  as  well  as  upon  the 
connectivity  of  the  pore  space.  We  find  Bedford's  [9]  argument  that,  in  some  local 
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sense,  each  constituent  will  obey  the  constitutive  relations  for  that  constituent  alone 
quite  appealing.  To  9et  up  macroscopic  constitutive  relations,  one  approach  would  be 
to  postulate  these  relations  directly  and  the  other  would  be  to  relate  macroscopic 
behavior  to  intrinsic  properties  of  the  constituents.  Apparently,  the  thermodynamic 
restrictions  on  entropy  production  as  well  as  the  notion  of  energy  density  are  applica¬ 
ble  to  each  constituent.  Bowen's  [24]  comments  regarding  mixtures  of  fluids  need  fur¬ 
ther  careful  investigation. 

Morland  [103-105]  did  not  assume  the  existence  of  an  energy  function  for  the 
mixture  but  still  admitted  interdependence  in  constitutive  relationships.  This  implies 
admitting  constitutive  relationships  with  possibly  nonsymmetric  coupling  effects  of  the 
type  proposed  earlier  by  Schiffman  [143]. 

Some  theories  involve  the  volume  fraction  of  the  constituents  in  addition  to 
the  usual  variables.  Morland  pointed  out  that  the  partial  density  variation  is  not  the 
density  variation  of  the  constituent  since  the  mixture  postulate  eliminates  reference 
to  the  actual  volume  occupied  by  each  constituent  in  an  immiscible  mixture.  A  correct 
theory  would  ensure  that  deformation  be  associated  with  a  set  of  particles  rather  than 
with  a  fixed  volume  in  space. 

For  a  statistically  isotropic  saturated  material,  Biot  expected  the  kinetic  energy  of 
the  saturated  soil  to  be  quadratic  in  the  velocities  of  the  fluid  and  the  soil  and  a 
coupling  term  was  included.  This  introduced  an  inertial  coupling  between  the  soil  and 
the  fluid.  It  is  difficult  to  assign  numerical  values  to  the  various  quantities  that  arise 
as  a  result  of  this  coupling.  Implementation  of  this  theory  in  a  finite  element  com¬ 
puter  program,  and  a  preliminary  parametric  study  to  investigate  the  efect  of  this 
coupling  on  the  dynamic  response,  indicated  that  the  effect  would  be  insignificant. 
However,  further  study  of  this  particular  feature  is  needed. 


It  appears  that  a  theory  of  dynamics  of  saturated  soils  should  use  a  convected 
coordinate  system  to  define  the  motion  of  the  soil  so  that  the  same  set  of  particles 
constitute  the  reference  volume.  The  flow  of  the  fluid  should  be  considered  as  relative 
to  this  reference  set  of  soil  particles.  Stress  would  be  defined  in  terms  of  these  con¬ 
vected  coordinates  and  the  balance  equations  would  then  be  written  for  the  reference 
set  of  particles.  This  would  represent  a  generalization  of  Gibson's  theory  of  nonlinear 
soil  consolidation  to  three-dimensions  and  also  to  include  inertia  effects.  For  slow  flow 
and  small  deformation,  certain  simplifying  assumptions  would  lead  to  Biot's  theory.  The 
soil  and  water  have  relative  motion  prior  to  liquefaction.  At  liquefaction,  the  relative 
velocity  reduces  to  zero  and  the  soil-water  mixture  would  move  as  a  single  fluid. 
Development  of  constitutive  equations  would  involve  the  thermodynamics  of  the  con¬ 
stituents  including  their  interaction  but  it  would  not  include  assigning  physical  mean¬ 
ing  to  a  "mixture  in  motion”.  Constitutive  and  inertial  couplings  might  exist.  How¬ 
ever,  symmetry  of  these  couplings  cannot  be  apriori  claimed  on  the  basis  of  the 
existence  of  energy  functions  for  the  mixture.  The  next  section  describes  such  a  theory 
developed  by  Hiremath  [75]. 


Section  DI 


A  DYNAMICAL  THEORY  OF  SATURATED  SOILS 

3.1  INTRODUCTION 

In  existing  theories  of  mixtures  the  multicomponent  mixture  has  been  regarded  as 
a  set  of  superposed  continua  in  motion.  The  mixture,  at  any  instant  of  time,  has  been 
defined  as  a  set  of  particles  constructed  by  superposition  of  constituent  particles.  In 
reviewing  theories  of  mixtures,  including  their  possible  relationship  with  mechanics  of 
saturated  soils  and  liquefaction  phenomena,  an  important  finding  was  that  the  notion 
of  the  mixture  as  a  continuum  in  motion  is  inadmissible  except  in  the  case  of  no  rela¬ 
tive  motion  between  the  constituents.  Liquefaction  is  primarily  caused  by  the  relative 
motion  of  soil  and  water  and,  therefore,  a  correct  theory  of  liquefaction  cannot  be 
derived  from  the  assumption  of  the  saturated  soil  being  a  mixture  in  motion  as  a  con¬ 
tinuum. 

Because  the  mixture  cannot  be  viewed  as  a  continuum  in  motion,  it  appears  inap¬ 
propriate  to  define  energy  functions  on  the  "mixture”  consisiting  of  centers  of  mass. 
This  implies  that  in  setting  up  constitutive  relationships  for  the  mixture,  one  cannot 
invoke  the  existence  of  an  energy  function. 

Some  investigators  (e.g,  [46,76,860  considering  the  special  problem  of  flow 
through  deformable  porous  solids,  have  attempted  to  write  the  balance  equations  in 
terms  of  relative  motion  and  porosity,  which  is  essentially  a  measure  of  relative  defor¬ 
mation.  It  would  appear  that  a  theory  based  on  balance  equations  written  for  a  refer¬ 
ence  set  of  particles  of  the  porous  solid  would  be  the  most  appropriate  for  this  case. 
Gibson  [55]  developed  such  a  theory  for  the  quasi-static  problem. 
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Hiremath  developed  a  [75]  theory  of  dynamics  of  saturated  soils  based  on  use  of 
convected  coordinate  system  to  describe  the  motion  of  a  fixed  set  of  soil  particles  in  a 
reference  volume  and  regarded  the  flow  of  water  to  be  relative  to  this  volume  of  the 
soil  (Item  4.5  in  Appendix  B).  This  theory  may  be  regarded  as  an  extension  of  the 
concepts  presented  by  Gibson  [55]  for  the  case  of  one-dimensional  quasi-static  deforma¬ 
tion  of  soils,  to  three  dimensions  and  to  include  inertia.  In  the  remainder  of  the  sec¬ 
tion  we  summarize  the  salient  features  of  this  theory. 
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33  BALANCE  LAWS 


33.1  Mass  Balance  of  the  Solid 

If  p^l)  and  pu)  denote  the  mass  densities  in  the  configurations  C,  and  C  respectively,  of 
a  volume  element  containing  the  same  set  of  solid  particles, 

JO  -  JO  JO*  (i 

P0  no  Po  (33) 


p(1,-n(,,p(1>  (3.4) 

in  which  n^1}  and  n(,)  are  the  solid  volume  fractions  in  the  initial  and  the  current  con¬ 
figurations,  respectively,  and  a  superposed  asterisk  denotes  an  intrinsic  quantity.  This 
leads  to  an  equation  of  mass  balance  in  the  form, 

fp^^dV,-  /V'V’dV  (3.5) 


V#  and  V  denote  volume  of  the  same  set  of  particles  in  Q,  and  C,  respectively.  Using 


(A.126), 


0<n1)*  n*11  —  s/G  p<1)’n(1)]  dV  =0 


The  point  form  of  this  equation  is 

JO*  JO  _  /fJ*>*  (O 
Pn  *  vGP  n 


If  the  solid  is  incompressible,  p<01>  =  p(1>*  and 
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3.3.2  Mass  Balance  of  the  Fluid 

The  motion  of  the  fluid  is  relative  to  the  solid  and  the  fluid  itself  does  not  have  a 
reference  state.  Mass  continuity  of  the  fluid  constituent  is  described  for  a  solid  vol¬ 
ume  in  the  current  configuration. 

Consider  a  rectangular  parallelepiped  of  the  solid  at  point  P,  in  the  reference  con¬ 
figuration  C*  which  corresponds  to  a  skew  parallelopiped  at  point  P  in  the  current 
configuration  C  The  parallelopiped  in  C,  is  made  up  of  surfaces  x,  -  constant  and 
x,  +  dx,  -  constant  (Figure  3),  and  during  motion,  encloses  the  same  solid  particles. 


Recalling  (A.110)  and  (A.112),  the  velocities  of  the  solid  and  the  fluid  in  terms 
of  the  base  vectors  e,  and  G,  (or  G1)  are 

t(1)  =  v(1)"G  =  v(,)Gm  =  u(1)e  (3.9) 

n  is  in  od 

and 


(3.10) 


The  components  u^1  (k  =1,2)  are  quantities  associated  with  the  reference  state.  The 
face  of  the  parallelopiped  formed  by  the  sides  dx3  and  dx5  in  the  reference  state 
becomes  an  area  formed  by  the  vectors  G,  dx5  and  G,  dx5  in  the  deformed  state.  Thus, 
the  mass  flux  per  unit  time  entering  this  face  is 

n'V”V|il-u(|,W<U!  (3.11) 

The  mass  flux  leaving  the  opposite  face  is 

a(V”Vl”  -  u'/W'dx1  +  -t.  {n'W,8  -  i'1  W)  <to2  dxJ  (3.12) 

The  net  rate  of  gain  including  the  fluid  flow  in  all  the  three-directions  is, 

(3.13) 

8*' 


The  rate  of  increase  of  fluid  within  the  deformed  parailolepiped  is 
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Ab(V*V5dx1dx2d*3] 

d* 


(3.14) 


Adding  (3.13)  and  (3.14)  and  cancelling  dx1  dx3  dxs,  the  equation  of  fluid  mass  continu¬ 
ity  is 


4-ln‘V8-  JSi  +  J- W"p<8,ti“-i|1>l)-0 


(3.15) 


dx 


3-3-3  Balance  of  Momentum  of  the  Fluid  Phase 

Consider  the  elementary  parallelopiped  that  is  rectangular  in  the  reference  config¬ 
uration  and  is  transformed  into  a  skew  parallelopiped  in  the  current  configuration  (Fig¬ 
ure  4). 


Let  —  t^J  denote  the  stress  vector  of  the  fluid  phase  acting  normal  to  the  surface 
formed  by  the  vectors  G2dxJ  and  G,  dx5.  The  net  force  due  to  this  stress  across  the 
surface  is,  noting  (A.132)  and  (A.142); 

- 1(2)  VgG11  dx2 dx3  =  ¥2) dx2 dx1  (3.1 6) 

The  internal  forces  on  the  six  faces  of  the  parallelopiped  then  are; 

-¥2)dx2dx\  TC2)dx2dx3  +  -^rT(J2>dx1dx2cli3 

dx 

-T^dx'dx3,  ¥2) dx 1  dx3  +  -&j  ¥22) dx 1  dx2  dx3  (3.17) 

dx 

—  T^2)dx!  dx2,  T(32)dx,dx2  +  -^-T<32)dxldx2dx3 

dx 

Similarly  the  net  body  forces,  inertial  forces  and  viscous  coupling  forces,  respectively, 
acting  over  the  deformed  volume  VG  dx1  dxJ  dx5  are; 

[p(2V2),  p(2V22),  p(2V32)]  VGdx1  dx2dx3  (3.18) 

[-p(2¥2),  -p,2¥22),  -p(2¥32>]  v/Gdx‘dx2dx3 


(3.19) 


T*,2,dx‘  dx2  +  T^'dx1  dx2  dx3 

a*3 
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{-D(t<12)-y(11)),  —  D (Vj2)  —  Yj1^,  -D(y(32)-y(31))}  •s/Gdx1dx2dx3  (3 JO) 

Summing  up  the  forces  and  setting  the  total  equal  to  zero  for  equilibrium, 

-3-  if }  +  Jg  p(2) F*2)  -  y/G  P(2)  f*2)  -  >/G  D [y(2)  -  y(,)]  =  0  (3.21) 

ex1  1 

Upon  use  of  (A.141),  and  rearranging 

_L_  JL [7g  r(2)ijGj  +  p(2V2)  =  p(2) f^2)  +  D(y(2)-  y(1))  (3J2) 

>/G  3x‘  r 

Recalling  (A.110),  (A.111),  (A.112)  and  (A.113)  along  with  (A.143)  gives 

t<81\  +  p<8  F181  =  p<!l  f*211  +  D  [v,8i  -  vl0)l  (M3) 


Alternatively,  using  (A.114)  one  can  refer  the  quantities  in  (3.22)  to  the  reference  state 
C„.  and  write 


-i-  [  s/'G  t(2*j  z  el  +  JG  p(2)  ^2)  e  =  n/G  p<2)  uC2>  e+  yGD  (ii(2)  -  u'*)  e  (3.24) 

i1,  m.jm  mm  “  m  mmm 

OX 

which  gives,  by  (A.150)  and  (A.151), 


MS*#?-  ^„(2,«<8+  VBD(i«-i«’) 


(3-25) 


or 

[S(2)i]  +  VGp(2)^2>=  v/Gp<2)u<2>+  VGDCu^-u^  (3J6) 

m  ,  i  m  in  m  o 

For  isotropic  fluid, 

S(2)'  =  yS‘  (3J7) 

n  n 

Then  (3.26)  gives, 

7T  +  yG  p(2)  ^2)  =  v/G  p(2)  u<2)  +  VG  D  (u(2)  -  d(,))  (3.28) 

,  Ill  10  uj  ID  ID 
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3J.4  Balance  of  Momentum  of  the  Fluid-Saturated  Solid 

Considering  the  reference  set  of  solid  particles  in  the  current  configuration  of  the 
deformed  parallelopiped  Figure  5,  let  —  t,  denote  the  total  stress  vector  acting  on  the 

strained  body  per  unit  area  formed  by  the  vectors  G3  dx2  and  G,dx’.  The  net  force 
across  this  surface  using  (A.142)  is, 

—  Tj  dx2  dx3  =  -t,  V  GG11  dx2  dx3  (3.29) 

The  other  qtiantities  viz.  T}  and  T,  are  defined  likewise  by  cyclic  permutation  and  are 
shown  in  Figure  5.  The  forces  on  the  six  bounding  surfaces  are: 

—  T,  dx2dx3 ,  T1dx2dx3  +  -3TT1dx‘dx2dx3 

3* 

—  T2dx*  dx3,  T2dx1dx3  +  -^-T2dx1dx2dx3  (3.30) 

3* 

—  T  dx1  dx2 ,  T1dx1dx2  +  -$7T,dx‘dx2dx3 

3x 

For  =  F(J)  =  F  the  body  forces  and  inertial  forces  acting  over  the  deformed  volume 
VCJdx1  dx3  dx1  can  be  expressed  as, 

[pF,t  pF2,  pFj]  -/Gdx1  dx2dx3  (3.31) 

and 

[-{p<,lf‘|,l  +  p<^,f‘1^,),  -V,V“  +  (>uV1a).  -fp'X'  +  p^lVadx'di'di5  (3.32) 
Summing  up  the  forces,  we  have, 

-$-T.+  VGpF=  jG[p{1)  +  p(2)  ?2}]  (3.33) 

3* 

Using  (A.141)  in  (3.33) 

_1_  -4- [>/g riJGj  +  pF  =  p(1) f*1*  +  p<2) f*2>  (3.34) 

VG  3x' 

Recalling  (A.62),  and  (A.110)  through  (A.113), 

r\  +  p  FJ  =  p(,Wl)i  +  p(2)  f U>i  (3.35) 
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Referring  quantities  in  (3.34)  to  the  reference  state  C0  and  using  (A.150)  and  (A.151), 
two  alternative  forms  of  (3-34)  are, 

[sij  z^j) ,  +  \/G  p  P  m  =  n/G  p(1) ««  +  ■JG  p<2)  u(m2)  (336) 

and 

[S'j,  +  y5pf.=  v/GpIU»l”  +  yOp(!,'i“  (337) 


3,1,5  Balance  of  Angular  Momentum  of  the  Fluid  Saturated  Solid 

In  absence  of  body  couples,  the  angular  moments  of  these  forces  with  respect  to 
the  deformed  axes  along  G,  are  given  by, 

-(TjXRWdx3.  (T1xR)dx2dx3  +  -^-(T1xR)dx‘dx2dx3 

3* 

-(T  xR)dx‘dx3.  (T, x R )dx1  dx3  +  -Ar (T, x R)dx‘ dx2dx3  (338) 

2  0* 

—  (T,  x  R)dx*  dx2,  (T3  x  R)dx*  dx2  +  (T3  x  R)dx*  dx2dx3 

Qx 

For  body  forces  and  inertial  forces  the  angular  moments  are,  respectively, 

IpF,xR,  pF2xR,  pFjXR]  >/Gdx’dx2dx3  (339) 

and 

[-V'V.'V’V,8).  -tp'X'+r'X'l- 

XR  VG  dx1  dx2  dx3  (3.40) 

Summing  up  the  moments, 

-$_(T  X  R)  +  (p  F  x  R)  >/G-{p(1)f(1)  +  p(2)f(2>}  X  R  JG  =  0  (3.41) 

3*'  ' 

or,  by  (A3  8)  and  rearrangement  of  terms 

[iH  +  jGpT-  jG{p(l}fll)  +  p(2)S2)]x*  +  TixG.i  =  0 

3* 

Noting  (333)  and  (3.42) 


(3.42) 
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T,  X  G,  =  0  (3.43) 

Using  (A.141),  this  gives 

>/G  r,j  G.  X  G(  =  0  (3.44) 

or,  equivalently 

€ljkriJ  =  0  (3.45) 

which  is  same  as 

r (3.46) 
Further,  in  view  of  (A.103), 

s'j  =  s^*  (3.47) 

The  bulk  stress  tensors,  r*J  and  s*J,  are  symmetric.  This  is  in  line  with  Green’s  [64] 
assertion  that  the  partial  stresses  need  not  be  symmetric  but  the  total  (bulk)  stress  is 
symmetric. 

3.4  SOME  SPECIALIZATIONS 

3-4.1  Specialization  to  One-Dimensional  Problem 

The  direction  z,  is  referred  to  simply  as  x  and  the  associated  quantities  are  denoted  by 
a  super-  or  sub-script  x. 


a).  {Cinema tical  Quantities 
From  (A.86)  and  (A.87), 


z  =z  (x,  r) 

(3.48) 

and 

x  (r  =  0)  =  a 

(3.49) 

For  the  motion  to  be  possible,  from  (A_5), 
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I-&I  >  0 

3* 


(3.50) 


The  position  vectors  of  points  P0  in  Q  and  P  in  C  are,  respectively,  ((A.77)  and  (A.78)) 

(331) 

(332) 


r=ie=ae 
R  =  ze 

The  displacement  vector,  n,  is  (A.79) 
u  =  R  —  r  =u  e 
u  =  z  — x  =  z—  a 

The  vectors,  G„,  Gx  and  metric  tensors  G„,  G“  may  be  defined  for  the  system  x 
in  the  body  at  time  t.  From  (A.97)  through  (A.99),  we  have  in  one-dimension, 

(335) 


(333) 

(334) 


G  =  R  =  £Le  =  (l  +  $-)e 
1  •*  3*  8* 


Gx  =  ^e 
8z 

o  =o,.G  =^^.=a  +  *-)2  =  l+2*.+(f)2 

“  *  *  3x  3*  3*  8*  o* 

G“  =  &-  &- 

dz  3z 

The  line  elements  are  given  by  (A.lOl)  and  (A.102), 


ds  =  G  dx  dx 


dSg  =  dx  dx 

The  strain  tensor,  using  (A.  100)  is, 

2  ^  **  3x  2  3x  0x 

The  changes  in  volume  are  obtained  from  (A.123)  and  (A.126)  as, 
dV0  =  dx 
dV=7GdV 


(336) 

(337) 

(338) 

(339) 

(3.60) 

(3.61) 

(3.62) 

(3.63) 


where. 
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G  =  IGI  =  1(1 +$-)2 1=  (1  +  ^-)2  =  (-&-)2 
“  Qx  Qx 

(3.64) 

(3.63)  and  (3.64)  give 

dV  =  (l+-3jj-)dx 

8* 

(3.65) 

The  expressions  for  velocity  and  acceleration  derived 

from  (A.110)  through  (A.123)  for 

one-dimension  are. 

v  *  u  e  =  v1  Gx  =  vx  Gx 

(3.66) 

f  =  ue  =  f*GJt  =  fxG* 

(3.67) 

b).  Mass  continuity  of  the  solid. 
From  (3.7) 


Ji>*  <U_  /r  Jll'Jl) 

P0  no  ='/°P  n 

(3.68) 

Using  (3.64),  mass  continuity  in  one-dimension 

(i)*  ( i)  _  to.  Air  At) 

0  n_  =  d  n 

00  ax 

(3.69) 

This  expression  is  the  same  as  Gibson's  [55] 

c).  Mass  continuity  of  the  fluid. 

From  (3.15) 

A  [n(2) p(2)*  To)  +  -3-  {n(2) p(2)> [vi(2)  -  d(,)] }  =  0 

at  a* 

(3.70) 

which,  upon  use  of  (3.64),  is 

A  [*<» p<2>*  JSL  {n(2)  p(2)*  [u(2)  -  i(  1 °] }  =  0 

at  a*  ax 

(3.71) 

This  is  the  same  as  Gibson's  equation  [55]. 


Momentum  balance  of  the  fluid. 


(3-28),  for  one-dimensional  analysis,  using  (3.64)  is; 

&L  +  4L  f+f*  =  u(2)  +  &.JL  n(2)  (u(2)  — u(I>)  (3.72) 

Jl  jx  *  $X  0x  ft 

If  the  inertia  term  is  neglected,  Gibson's  [55]  equation  for  quasi-static  analysis  is  recov¬ 
ered  viz. 

i2L  +  $L  pW&»  =  &.£  n(2)(u<2)  — u(1>)  (3.73) 

$x  flx  1  §x  « 

e).  Momentum  balance  of  the  fluid-saturated  solid. 


Recalling  (3.64),  for  one-dimension,  (3J7)  is 


&pp  =^.p(,)u(14-^-p<2)u(2> 

6*  x  0*  d* 


(3.74) 


Ignoring  the  inertial  term,  we  recover  Gibson's  [55l  equation  of  motion  for  the  bulk 


viz^ 


+  p  f  =  0 
1  0* 


(3.75) 


3.4.2  Small  deformation  Theory 

Biot's  [17,19]  equations  for  small  deformation  theory  are  embedded  in  the  general 
theory  presented  in  this  work  as  a  specialization.  Assuming  small  strain,  explicit  forms 
of  continuity  equations  are  not  required  as  the  changes  in  density  are  small  All 
quantities  are  referred  to  the  initial  state  with  rectangular  cartesian  system  as  a  frame 
of  reference.  In  that  case,  the  distinction  between  the  contra  variant  and  covariant  com¬ 
ponents  dispppears.  The  kinematical  relations  in  (3.1)  and  (3.2)  reduce  to 
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The  momentum  balance  equations  in  terms  of  the  bulk  stress  ty  and  the  partial  pres¬ 
sure  ir  are  obtained  from  (3-37)  and  (3.26),  respectively,  as 

ty  j  +  pF.  =  p(1)uj1)  +  p(2)u|2)  (3.77) 

tt.  +  p(2)  F.  =  p(2)Uj2)  +  D[Uj2)  -  ujl>]  (3.78) 

which  is  often  also  written  in  the  form; 

V, + pnr  f, = P<M"iu>  t£.a  -  d'11]  0.79) 

,i  I  i  K  1  1 

The  above  equations  are  the  same  as  in  one  formulaion  of  Biot's  theory.  Subtracting 
(3.78)  from  (3.77),  an  equilibrium  equation  in  terms  of  the  partial  solid  stress  is 
obtained  viz. 


t(1)  +  p(1) F.  =  -  *-  n(2)  n<2)[u<2> -  ii(1>] 

IJ.J  I  1  K  11 


(3.80) 


Comparing  with  Biot's  £17,18]  equations,  (3.79)  and  (3.80)  do  not  have  the  mass 
coupling  terms. 


3-5  CONSTITUTIVE  RELATIONS 

The  issue  of  defining  mechanical  quantities  for  which  constitutive  relationships 
need  to  be  defined  has  been  discussed  in  Section  II  and  in  the  technical  report  listed  as 
item  12  in  Appendix  B.  The  dynamical  theory  summarized  in  this  section  is  based  on 
studying  the  movement  of  a  connected  set  of  non-interpenetrating  particles.  In  this 
theory,  the  constitutive  equations  are  required  only  for  the  partial  stresses  and  the  dif¬ 
fusive  resistance.  Porosity  is  a  quantity  directly  related  to  the  deformation  of  the  set 
and  need  not  be  treated  as  an  additional  variable.  The  relative  movement  between  the 
pore-water  and  the  reference  set  of  soil  particles  would  apparently  be  the  principal 
kinematical  variable  related  to  the  interaction  force.  For  linear  theory  this  relationship 
would  reduce  to  d'Arcy's  rule.  The  stresses  in  the  reference  set  of  particles  must  be 
described  in  convected  coordinates  Just  as  the  deformation  and  defromation  rates  are. 


For  the  case  of  rate-independent  materials,  a  theory  for  a  single  elastic-plastic  material 
was  presented  by  Ayoub  [8J.  He  used  the  Cauchy  stress  as  the  mechanical  variable. 
He  showed  that,  for  large  deformations,  the  difference  between  the  conventional 
description  of  stress/strain  relaions  using  quantities  referred  to  the  original  configuration 
and  the  correct  description  proposed  by  him  would  be  quite  significant.  The  procedures 
and  descriptions  suggested  by  Ayoub  can  be  easily  generalized  to  admit  possible 
coupling  between  the  partial  soil  stresses  and  the  fluid  pressures.  There  is  apparent 
need  for  the  development  of  data  on  the  behavior  of  saturated  soils  under  large  defor¬ 
mations  to  define  the  nature  of  the  constitutive  relations.  Morland’s  [103]  proposal  that 
the  constitutive  equations  for  effective  stresses  in  the  porous  material  be  assumed  to 
have  the  same  form  as  that  for  the  intrinsic  material  appears  to  be  attractive  but 
needs  verification.  The  distinction  between  the  partial  stress  and  the  effective  stress 
would  allow  for  the  possible  coupling  between  the  constitutive  equations  for  the  soil 
and  the  pore-water. 


Section  IV 


SOLUTION  PROCEDURES 


4.1  INTRODUCTION 

The  solution  procedures  for  the  initial  value  problem  of  dynamic  response  of  soil 
imew  can  be  classed  into  the  following  groups, 

1.  Exact  Solutions 

2.  Semi-Discrete  Solution  Procedures 

3.  Finite  Element  Solutions 

Exact  solutions  were  developed  for  the  linearized  version  neglecting  mass  and  con¬ 
stitutive  couplings  and  assuming  that  the  water  was  completely  free  to  move  relative 
to  the  soiL  This  essentially  implied  a  specialization  to  Biot’s  theory.  The  exact  solu¬ 
tions  described  in  items  1.14,  1.15,  2.7,  2.8,  3.5  in  Appendix  B  and  the  semi-discrete 
methods  described  in  items  1.8  and  42  of  the  same  Appendix  were  based  on  this  theo¬ 
ry.  For  the  purpose  of  numerical  solution  (items  13,  1.6,  1.10,  1.12,  1.13,  1.17,  2.3,  2.4, 
3,3,  3.4,  4.i  in  Appendix  B),  nonlinearity  and  couplings  could  be  accomodated  to  a  cer¬ 
tain  extent.  In  this  section  we  describe  briefly  the  results  of  the  research  under  each 
of  the  three  headings. 
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4.2  EXACT  SOLUTIONS 

43.1  Introduction 

Exact  solutions  to  Biot's  equations  of  dynamics  of  fluid-saturated  porous  media 
were  obtained  by  Biot  [17-19].  Later  Deresiewicz  [40],  Chakraborty  [33l  and  Garg  [50] 
obtained  solutions  for  various  boundary  conditions.  In  the  present  research  the  work 
was  extended  and  computer  codes  for  numerical  solution  were  tested  against  these  exact 
solutions.  Items  1.14,  1.17,  2.7,  2.8,  and  3.5  in  Appendix  B  contain  details  of  this 
development.  The  specific  items  of  research  included  the  following: 

a) .  Garg's  fundamental  solution  for  the  problem  of  one-dimensional  wave  propaga¬ 
tion  in  fluid-saturated  media  was  integrated  to  develop  solutions  for  several  cases  of 
surface  loading  of  a  saturated  soil  column  of  infinite  extent.  [Items  1.14  and  2.7, 
Appendix  Bj. 

b) .  In  order  to  obtain  a  solution  to  the  problems  of  "strong  coupling”  "weak 
coupling”  Garg  had  made  certain  assumptions.  These  assumptions  were  carefully  inves¬ 
tigated.  [Item  3.5,  Appendix  B]. 

c) .  Solutions  to  Biot's  equations  of  wave  propagation  involving  sudden  changes  in 
excitation  were  developed  by  separating  the  propagation  of  the  singularity  from  the 
diffusive  process.  These  solutions  were  extended  to  some  two-dimensional  cases.  [Item 
1.15  and  2.8,  Appendix  B). 

In  the  following  paragraphs  we  summarize  some  of  the  results  of  these  investiga¬ 


tions. 


79 


4 


4JL2  Garg's  Solution 

Figure  6(b)  shows  the  load  for  which  Garg  obtained  an  exact  solution  and  the 
four  load  cases  for  which  additional  solutions  were  obtained  as  part  of  the  present 

research.  Garg  [13]  wrote  Biot's  equations,  for  the  one-dimensional  problem,  without 

inertial  mass  coupling  in  the  form; 

p(1)u  =  au  u  +  cU  u  —  D(u  —  U)  (4.1) 

p(2)U  =  cu  u  +  bU  n  +  D(ii-U)  (4.2) 

where  u,  U  are  the  displacements  of  the  solid  and  fluid  respectively.  Garg  [50] 
assumed  the  displacements  of  the  constituents  to  be  specified  on  the  end  x  —  0  as 


u(0,t)  =  f(t) 
U(0,t)  =  g(t) 


(4J) 


Only  the  conditions  at  x  =  0  were  needed  as  the  column  was  assumed  to  be  of  infinite 
extent.  The  initial  conditions  were; 
u  (x ,  0)  *  u0(x) 

(4.4) 

u  (x ,  ,0)  =  u0(x) 


U(x.O)  =  U0(x) 


CJ(x.O)  =  UQ(x) 

In  order  to  solve  the  wave  equations,  (4.1)  and  (4.2)  were  differentiated  with 
respect  to  the  time  variable.  For  homogeneous  initial  conditions,  the  boundary  conditions 
on  the  velocities  of  the  constituents,  assuming  they  move  together  at  this  point,  were 

u(0,t)  =  U(0,t)  =  <6(0  (4.5) 

Applying  the  Laplace  transform  to  the  time  derivatives  of  (4.1)  and  (4.2),  and  denoting 
the  velocities  of  the  solid  and  the  fluid  by  v,  V, 


(l)pzv  =  a  -3—^-  +  c  3  +  Dp(V-v) 


d* 


* 


J 


(4.6) 


81 


p(V  v  =  ciii  +  b^.  -  Dpcv-v) 

ax2  ax2 


where 


L[u(x,t),U(x,t)3  =  [v(x,p),V(x,p)] 


indicates  Laplace  transformation,  and  p  is  the  transform  parameter.  Assuming  the  solu¬ 
tion  to  have  the  form 


v  _  A 
V  ~  B 


exp(  — yx) 


the  characteristic  equation  is: 


>'p(cy-p2)-(cty2-p2)(civ2-p2)  =  0 


(4.10) 


where: 


02  _  a  +  b  +  2c 
C°  P 


2C2  =  c\  +  C*  ±  [(c;  -  C +  4Cj2C^] 
cj  =  a/p(0 

c2  =  b/p(2) 


<2\2  ,  ~2  il/2 


(4.11) 


_  /  JO 

Cl7  -  C/p 


C',  =  c/p(: 


JO  (2) 

P  P 


C0  is  the  wave  velocity  when  the  saturated  medium  acts  as  a  single  material  and  C± 


are  the  wave  velocities  when  no  viscous  coupling  exists.  (4.10)  has  four  roots,  viz_ 


y\  2  =  M,(p)  ±  [m;(p)-m2(p)]‘ 


(4.12) 


where 


2M,(p)  =  (l/C±  +  1/C2)p2  +  vpC2/(cfc2) 
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M2(p) 


p3  (  p  +  v  ) 

c*ci 


Thus,  the  general  solution  can  be  written  as 


v(x,p)  =  A,  exp(— ytx)  +  A2  exp(-y2x) 
V(x,p)  =  Bj  exp(— yjx)  +  B2  exp(—y2x) 


(4.13) 


Here  A1,  Aj,  B,  and  B,  are  functions  of  p  and  must  be  determined  by  the  following 
boundary  conditions  and  compatibility  equations. 

A,  +  A2  =  ?(p) 

B,  +  B2  =  ?(p) 

[p2  -  C\y\  +  Dp/p(1)]A,  =  [C\2y\  +  Dp/p<1)]B,  (4.14) 

[p2  -  C\y22  +  Dp/pU)]A2  =  [C\2y22  +  Dp/p(l)]B2 


Hence, 

Aj  =  ?(p)(l-S2)/(S1-Sp 

a2  =  ^(P)(sl-i)/(s1-s2) 


(4.15) 


where 


s,  =  (P2-cjy;  +  Q)/(c;2V;  +  Q) 
s2  =  (p2~C|yj  +  Q)/(Cj2y2  +  Q) 


A  general  solution  based  on  inverse  transformation  of  (4.13)  is  not  available.  Two  spe¬ 
cial  cases  were  solved  by  Garg  [50]  For  relatively  small  value  of  D,  Garg  approximated 
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7)2  =  ( C2Q  -  C2J/(.C2+  -  C2_ )  (4.19) 

* 

As  D  approaches  zero,  7^,  vanish,  and  the  expressions  for  amplitudes  in  (4.15) 
reduce  to; 


A,  =  A*$Kp) 

A2  =  A*?(p)  (4.20) 

Bj  =  b;?(P) 

B2  =  B*$(p) 

where 

a’  *  (c{  -  cf.  +  c212)/(cl  -  cl) 

b]  =  (c*  -  c2  +  c2n)/(cl  -  ci) 

A*  =  1  -  A*  ,  b'2  =  1  -  b‘ 

Garg  termed  this  special  case  "weak"  coupling.  In  evaluating  A*,  B|,  he  set  D  =  0,  i^, 

no  coupling  to  get  7,  =  p/ C*  and  y,  =  p / C_ .  This  assumption  to  avoid  dependence  of 
the  amplitudes  on  the  transform  parameter  p  made  inversion  of  the  solution  possible, 
but  is  inconsistent  with  the  approximation  to  get  equations  (4.15)  and  (4.16).  Substitut¬ 
ing  these  results  into  (4.13),  the  transformed  solution  for  weak  coupling  was  written 


as; 
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$  (p)  exp  [  -  (x/C+)  F ,  ] 


+ 


$  (p)  exp  [  —  (x/C_)  F2  ] 


(4J21) 


where 

Fj  =  [Cp  +  t,,)2-^]1'2 

F2  =  Up  +  r,2)2-^l,,! 

Inversion  of  (4-21)  gave 

.‘I  [  exp  (-  rjjt)  8  (t  -  -p-)  +  exp  (-  rjjt)  ft(t)  H(t  -  ]  *  <(>  (t) 


u(x,t) 

U(x.t) 

B, 


[exp  (-  r)2t)  8(t  -  exp  (-  r)2t)  f  2(t)  H(t  -  ^-)  ]  *  <j>  (t) 


where 

_  ijytW/cm 

1  (t2  — x2/C2  )1/2 

_  I.t^-xW2] 

2  ( t2  —  x2/C2  )1/2 


(4.22) 


Here,  H(t)  is  the  Heaviside  step  function,  S  (t)  is  the  Dirac  delta  and  I,  is  the  modified 
Bessel  function  of  first  kind  of  order  one.  The  symbol  *  denotes  convolution  product. 
In  evaluating  the  amplitudes  no  viscous  coupling  was  assumed  while  effect  of  viscous 
coupling  approximated  to  the  first  order  was  retained  in  the  exponentially  decaying 
terms.  If,  for  consistency,  we  set  D  =  0  in  the  exponential  decay  terms,  the  trans¬ 
formed  solution  (4.21)  would  reduce  to: 


v(p)| 

K 

V(p)| 

k 

?(p)  exp[  — (x  p/C+)] 


+ 


?(p)  exp[  — (xp/C_)] 


(423) 


The  inverse  is: 


u(x.t)| 

A'. 

• 

8  (t  —  i/CJ  *  <>  (t)  + 

A’r 

• 

U  (x ,  t)j 

B. 

b2 

8  (t  —  x/C_)  *  (t) 


(4.24) 
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This  is  the  solution  for  no  viscous  coupling  in  which  the  solid  and  the  fluid  particles 
move  independently. 

As  D  goes  to  infinity,  the  characteristic  equation  (4.10)  yields  a  single  root: 
y  =  (4.25) 


which  corresponds  to  wave  propagation  with  speed  Q,  iA,  the  mixture  moves  as  a  sin¬ 
gle  material.  For  moderately  large  value  of  D,  Le,  "strong”  coupling,  Garg  [50]  wrote 
the  first  order  approximation  of  the  characteristic  root  in  1/D  as: 


y  =  JL[ _ l _ ] 

y  C0li  +  epj 


(4.26) 


where 


0  =  -  -McJ/cJ-utci/cJ-i] 


(4.27) 


The  expression  for  0  given  in  [50]  is  in  error.  (4.27)  is  the  corrected  form  given  by 
Garg  [51].  (4.26)  describes  the  motion  of  the  mixture  in  which  v  and  V  are  different 


in  order  1/D.  Hence,  Garg  [50]  assumed: 
v  =  V 


(4/28) 


Based  on  this  approximation,  the  transformed  solution  was  obtained,  using  equations 
(4.26)  and  (4.28)  as: 


7(p)  - 


(4/29) 


or,  equivalently: 


v (p)  =  ?(p)  exp[--ij]  exp[ - j~— 

Co0  Cn02(p  + 


p  + 1/9) 


(4.30) 


The  inverse  is: 


u(x,t)  =  exp[-  jr~a  ]  0(t)* 


+  8<0) 


(4.31) 
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4.23  Integration  of  Garg's  Solution 

Garg's  solution  (424),  for  "weak  coupling”,  was  integrated  to  obtain  explicit  solu¬ 
tions  for  four  different  velocity  boundary  conditions  shown  in  Figure  6.  These  solu¬ 
tions  are  listed  below.  Details  are  given  in  item  1.14.  of  Appendix  B. 

a).  Unit  Box  Function 


The  applied  velocity  boundary  condition  at  x  =  0  is 
0(t)  =  H(t)  -  H(t-t,)  (4.32) 

where  tj  is  the  time  at  which  the  excitation  is  reduced  to  zero.  Substitution  into 
(421)  gives 

{  exp  C  —  1  +  J  «p[“V]  f,(r)dr} 


x  H(t--i-)  -  H(t  — tj  —  -p-  ) 


B 


{exp[--^]  -  J'  exp  [  —  t)2t  ]  f  2(t)  dr  } 


x  H (t  —  )  -  H(t  — tj  —  •£-  ) 


(4.33) 


b).  Sine  Function 


Assuming  that  the  velocity  specified  on  the  boundary  is  harmonic;  Le^ 
^(t)  =  sin  (ait) 

where  a i  is  the  frequency,  the  corresponding  solution  is  given  by 

{sinb^t  —  -5-)  exp[—  ~ — [sin  bit  /*  fj(r)  cosbit  dr 
c+  C+  Jfc 


(4.34) 


—  COS  bit 


f  ,(r)  sin  bit  dr )  }  H  (t  —  x/C+  ) 


i 


• 

i 

i 

j 

I 

i 


J 
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+  2  { sin <u(t  —  ■— )  exp [  —  [sin  cut  f 

B2  -  -  C-  j7C 


f2(r)  coscut  dr 


—  coscut  I  f,(r)  sincut  dr] }  H(t  — x/C_) 


(4.35) 


c).  Ramp  Function 


Velocity  function  specified  as  acting  on  the  boundary  is 
<t>( t)  =  —  [H(t)  -  H(t  — tj)]  +  H(t-tj) 


The  solution  for  this  case  is: 


(4.36) 


+  f  exp[-T,1T]dT]H(t-x/C+) 

x7C+ 

+  [i-(t,-t  +  .p-)  exp[-^ii.] 

j  t 

+  f  f ,  (r)  ( 1  exp[  —  7),  r  ]  dr  ]  H (t  -  tj  -  x/C+)  } 

170 


2  1 


+  J  f2(r)~ —  exp[-772r]  dr]  H(t-x/C_) 

+  i-M, 

+  /*  f ? (t)  (1  — - — —  expl-rj,  r ]  dr]  H(t  — tj  —  x/C_ 

-  ,7c  li 


)  )  (4.37) 


88 


d).  Spike  Function 

The  excitation  in  this  case  may  be  expressed  as 
*(t)  =  —  [H(t)  -  H(t-t,)]  +  (  — -2)  [H(t-2t,)  -  H(t-t,)]  (4.38) 

The  solution  for  velocity  is: 

U1 

U| 

x  [H(t  — x/C+)  -  H(t  — tj  —  x/C+)]  +  [  J-(t--|--2t1)exp[-^] 


B, 


{[  —  (t-p-)  expt--^-]  +  f  f, (r)  ^  exp[— rjj t ]  dT ] 

W  W  W  dc.  i 


+ 


+  f  fj(r)  — — —  -  2)  exp[-7),  t]  dr] 

x  [H(t  — 2tj  —  x/C+)  -  H(t  — tj  —  x/C+)  ]  } 
{[i(t-±)«p(-^i]  +  ^  J  exp[— t)2t]  dr] 

x  [H(t-x/Cj  -  H(t  — tj  —  x/C_)]  +  [  J-(t-^--2t,)exp[-^-] 

+  /*f2(r)fcll-2)exp(-r)2r]dT] 

C-  it.  li 

x  (H(t-2t,-x/C_)  -  H(t  — tj  — x/C_)]  }  (4J9) 


4JL4  Evaluation  of  Garg's  Approximations 

To  obtain  exact  inverses  to  the  "weak”  and  the  "strong”  coupling  problems,  Garg 
[SO]  made  some  assumptions.  Primarily  these  amounted  tto  neglecting  the  Laplace  trans¬ 
form  parameter  p  in  comparison  with  D  or  its  reciprocal  depending  upon  whether  D 
took  on  very  large  (strong  coupling)  or  very  small  (weak  coupling)  values.  The  reason¬ 
ableness  of  these  assumptions  was  examined.  As  the  range  of  values  of  the  trans- 
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formed  parameter  p  extends  over  the  entire  positive  interval,  it  would  appear  improper 
to  compare  the  parameter  D,  or  its  inverse,  with  p  in  terms  of  order  of  magnitude. 
However,  if  0^  =  0,  =  C_,  the  solution  given  by  (431)  for  strong  coupling  is  correct. 
But,  for  that  case,  9  =  0  and  the  solution  reduces  to  the  one  for  the  case  D  -* «.  This 
special  case  requires  that  the  quantity 

(C2  -C^)2  +  4Cl22C2, 

should  vanish,  Le.  C,  =  C3  and  Cu  =  CJ1  =  0.  This  in  turn  requires  a/pu)  =  b/pa)  and 
c  =  0,  Le.  the  only  coupling  in  (4.6)  and  (4.7)  is  through  the  viscous  coupling.  Vanish¬ 
ing  of  c  implies  that  the  constitutive  relations  for  the  fluid  and  solid  partial  stresses 
are  uncoupled. 


For  the  case  of  weak  coupling  too,  Garg  used  a  linear  approximation  in  D  for 
roots  of  the  characteristic  equation;  but  for  determination  of  the  amplitudes,  D  was  set 
equal  to  zero.  If  the  same  linearization  is  used  for  the  amplitudes  as  well,  the  solution 
would  be 

Aj  =  -#p)J,/F,  A2  =  #p)J/F 

(4.40) 

B,  =  ^(p)L,/F,  B2  =  #p)L/F 

where 


J,  =  (C2 jA ,  +  D/Pg1’ )  —  (p  +  k) Q,  - 


C2  +  C2 


Jz  *  (Cj2AI  +  D/p“))-(p  +  K)Q,  -  "  D/Po* 

l,  *  P+(P+y)Qj +  (i-Sl£ii)D/p;,)-c;A1  ~(c;+c|px2 
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F  =  [C^p  +  Ccf  +  C^D/p^R,  +  C2l2R2 

R,  =  (CJ-Cj)/(C*Cf.) 

r2  =  [cjc^+ci)-2c^ci]/[c^cl(cj-ci)] 


X 


1.2 


Q. 


<a'p 

g^+c^) 

c2ci 


Q,  - 


^(cj+cjp 

dc2 


(4.41) 


The  amplitudes  are  dependent  on  p  in  a  complicated  fashion  and  an  analytical  inverse 
is  not  available.  Of  course,  for  v  -*  0  the  expressions  are  identical  to  those  in  [501.  To 
examine  the  applicability  of  Garg's  solution,  the  amplitudes  for  two  materials,  i.e,  the 
one  used  in  [50]  and  a  coarse  sand  with  mechanical  properties  listed  in  Table  1  were 
evaluated.  In  the  table,  K,  K,.  K,  are,  respectively,  the  bulk  moduli  of  the  porous  sol¬ 
id,  the  non  porous  solid  and  the  fluid,  p’  is  the  shear  modulus  of  the  non  porous  solid. 
(4.40)  gave  for  Garg's  material 

J,  «  —  (OJ238  p  +  101.45) 

J,  ■  -  (0.0044  p  +  1.455) 

L,  ■  -(0.1054  p  +  65.01) 

L,  »  (0.1139  p  +  34.98) 

F  «  (0.8906  p  4  5757.4) 
and  for  coarse  sand 

J,  «  -(54.966  p  +  148004.4) 

]}  =  -  (0.283  p  ♦  752.9) 

L,  =  -(1.4065  p  +  7386.2) 


4 


1 


•I 


92 


L,  -  (5X276  p  +  139863-6) 
F  -  (1497  p  ♦  57574) 


The  above  expressions  show  that  the  contribution  of  the  tens*  containing  p  is  negligi¬ 
ble  in  comparison  with  the  constant  terms  unless  p  takes  an  extremely  Urge  values. 
We  note  further  that  far  both  the  materials,  the  contribution  J,  of  the  second  wave  in 
the  solid  to  the  total  response  is  relatively  nail  For  the  coarse  sand  the  contribution 
L,  of  the  fust  wave  to  the  fluid  ■  also  relatively  small 


Derg's  approximate  solution  was  hardly  distinguishable  from  the  numerical  inverse 
o f  the  exact  transformed  solution-  It  seems  appropriate  to  conclude  that  Garg't  approxi¬ 
mate  solution  for  weak  coupling  is  acceptable  for  a  short  time  range  after  sudden 
application  of  uniform  velocity  at  the  end  of  the  column.  For  the  case  of  strong 
coupling,  it  appears  reasonable  to  set 

y  -  -?-(l-8)  (442) 

for  sufficiently  Largs  values  of  D.  (442)  can  be  rewritten  sa 

C^y1  -  p2  -  pW-26)  (443) 

and  results  in 


CW  -  p*[(^--l)  ♦  ^-(«a-28)) 

Co  Co 

Cy-p2  «  p^(— 7-1)  +  — ^  d2  -  25)) 

cj  cj 

Substituting  these  relations  into  (4.10), 


(4.44) 


(4.45) 


-2vpJ  *  P4(§-1)(%-1)  -  28(^(^.-l)  +  ^-(^--1)]  +  CXS2)  (4.46) 

/•<*  /^a  a  /~»<c 


C2  C2 

S>  ^0 


c  c 

^0  ^0 


•J 


where  (X81)  represents  terms  of  second  order  in  8.  Thus,  a  first  order  approximation  to 

8  is: 


where 


and 


J  ■  »  -  pA 


c*  c*  &  c* 

q  q  q  q 

For  the  above  result  to  match  Garg's  [SOI  it  ta 
lect  the  terms  which  contain  p.  (44*)  gave 

6  ■  131.6x10*  +  0J7?785p 
and  for  a  fine  sand 


(447) 


(448) 


to  assume  fi  ■  w.  ix,  to  neg- 


8  -  28541x10*  ♦  1.184847  p 

Garg'a  linearization  of  8  and  the  further  approximation  to  obtain  an  analytical  solution 
were  seen  to  yield  results  of  acceptable  accuracy  and  form  a  proper  basis  for  develop¬ 
ments  of  solutions  for  other  boundary  conditions. 


4JL5  Wavs  Propagation  in  a  Fluid-Saturated  Soil  Layer 

During  the  present  research,  the  coocepts  described  above  were  extended  to  a 
finite  uni  column  (elastic  sol  layer)  and  to  independently  specified  boundary  conditions 
for  the  two  constituents.  The  solution  procem  consisted  of  constructing  singular  fields 
which  incorporate  all  discontinuities  of  the  velocity  fields  and  their  first  and  second 
derivatives.  This  additive  decompaction  left  twice  continuoaly  differentiable  fields 
which  satisfy  coupled  hyperbolic  second  order  differentia]  equation*  with  continuous 
forcing  terms.  The  results  of  the  analysis  were  compered  with  numerical  inversion  of 
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the  i  transform  solutions  and  »i«o  with  numerical  solutions  obtained  by  using 

the  finite  element  method.  Details  of  this  research  are  documented  in  items  1.1S  and 
2-8  in  Appendix  B. 


4J  SEMI-DISCRETE  SOLUTION 

In  order  to  separate  the  error  contribution  of  the  spatial  and  the  temporal  approx¬ 
imations.  an  eigenfunction  approach  to  solution  of  spatially  discretized  equations  of 
mooon  was  developed.  This  was  expected  to  provide  a  benchmark  for  evaluation  of 
fully  discretized  time  dnmam  solution  procedures.  The  studies  showed  the  presence  of  a 
high  frequency  spurious  oscillatory  component  related  to  the  spatial  mesh  size.  With 
refinement  at  the  mesh,  the  lowest  eigenvalues  of  the  Laplace  transform  solution  were 
found  to  converge.  A  Ritz  vector  type  approach  in  which  the  base  vectors  are  related 
to  the  cxatatioo  could  possibly  improve  the  accuracy  at  the  solution. 


The  procedure  employed  consisted  of  s  finite  element  discretization  of  the  coupled 
aquations  of  motion  to  get  the  matrix  equations 


(449) 


Here  u.  w  are  the  soil  displacement  and  the  displacement  of  the  fluid  relative  to  the 
soil  respectively.  Cm  represents  solid  damping.  M,  K  and  C  are,  respectively,  the  mass. 


stiffness  and  damping  with  the  subscripts  at,  ff  and  if  indicating  the  respective  solid, 
fluid  and  coupling  components.  The  solid  damping  was  introduced  by  Ghabouan  [53] 
ss  a  linear  combination  of  the  stiff  ness  end  mass  matrices,  (Rayleigh  damping),  in  the 
following  fora; 

c„  -  .,(M.  -  f'M„)  +  -  «\,)  <«0) 

where  a,,  a,  are  constants.  (4.49)  can  be  expressed  compactly  as 
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[M]  {U>  +  [CHU}  +  IK]{U}  =  {R}  (451) 

where 


tad  [Ml  [Cl  [X]  tie,  respectively,  the  system  mie,  damping  and  tuff  acta  matrices. 
Dropping  the  square  brackets  for  convenience,  these  equations  after  Laplace  transforma¬ 
tion  becomes 

( s*M  ♦  sC  +  K  )q  «  F  (4.52) 

where 

q  ■  C  ■  C(s) 

(453) 

F  »  F(s)  »  R-f  sML'0  ♦  MU0  ♦  OJ0 


and  U, .  U,  are  the  initial  conditions  for  nodal  point  displacement  and  velocities. 
(452)  is  a  system  of  quadratic  equations.  To  Unesrue  the  system,  let 


q  *  sq  (454) 

Premultiplying  by  M,  and  rearranging, 

Mq  -  sMq  «  0  (455) 

Hence  (452)  can  be  written  as 

-(K  +  sOq-sMq  ■  — F  (456) 

Combining  (455)  and  (456): 

IM-MItl-H 

or. 

[A  -  sB]$  =  -f  (458) 


where: 


(4.58)  is  a  system  of  linear  equations.  A  and  B  an  symmetric  matrices.  The  pre¬ 
scribed  boundary  conditions  on  displacements,  velocities  etc.  wen  enforced  following 
Wilson's  method.  To  obtain  a  solution  to  these  aquations,  the  vector  4  was  expressed 
as  a  linear  combination  of  a  set  of  independent  vectors  Q v  n  “  1.2. _ m,  in  the  follow¬ 

ing  form; 

4“I>.  0. -10)U)  (4.59) 

when  a.  an  coefficients  and  Q.  is  the  n“  column  of  the  matrix  [Ql-  The  eigenvec¬ 
tors  of  the  problem 

[  A  -  sB  ]  y  -  0  ;  n  ■  \X—  m  (4j60) 

wen  taken  to  be  the  independent  vectors.  The  eigenvalues  a,  were  determined  as  the 
roots  of  the  polynomial  equation 

I  A  -  sBl  -  0  (4.61) 

(4J8)  and  (4J9)  give 

(  A  -  sB)Qa  -  -f  (4,62) 

Premultiplying  both  sidei  by  QT 

A  a  -  sB  a  ■  -QTF  (4.63) 

when  A*  *  QTAQ  and  B*  *  QTBQ  an  diagonal  matrices  due  to  orthogonality  of  Q  with 
respect  to  A  and  B.  The  nu  equation  is 

A\  ~  *  -Ql1"  ( '<*«> 


where 
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Q>Q. 


,Q>5. 


*.B. 


(4^65) 


Substitution  of  (4-65)  into  (464)  gave 

■X\  -  *>.  -  -0? 

Hence 

-c£r 

a  ■  - - - 

‘ 

Substitution  of  (467)  m  (459)  gave  the  solution  in  the  Laplace  transform  space; 


(466) 


(467) 


A  Qlr 

*-,(s-s.)B 


(468) 


This  series  was  inverted,  term  by  term,  to  get  the  required  results  u>  the  form  of 
nodal  displacements  and  velocities  as  functions  of  time.  These  functions  were  evaluat¬ 
ed  for  specified  values  of  the  time  variable  to  determine  the  solutions  as  well  as  sec¬ 
ondary  quantities  of  interest  eg,  the  stress  in  the  matehaL 


Items  1.8  and  42  in  Appendix  B  contain  details  of  the  procedures  as  well  as 
examples  for  validation  of  the  computer  codes. 


44  FINITE  ELEMENT  SOLUTIONS 

44.1  Variational  Formulation 

V* national  formulations  of  the  linearized  version  of  the  theory  and  its  extension 
to  include  material  nonlinearity  were  developed  to  construct  a  basis  for  alternate  finite 
element  approaches  to  the  problem.  Ghabouasi  [S3]  had  developed  a  variational  formula¬ 
tion  of  Biot's  theory  but  for  the  purpose  of  finite  element  analysis  he  used  the  Galer- 
lun  procedure.  Also,  this  variational  formulation  did  not  allow  for  the  boundary  con¬ 
ditions  properly  nor  did  it  allow  for  interelement  discontinuities  inherent  in  finite 
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element  method*. 


In  order  to  systematically  develop  variational  principle*,  the  equation*  of  motion 
were  written  in  *  form  that  they  would  constitute  a  «e  If -adjoint  system  in  an  appro¬ 
priate  linear  vector  apace.  Thi*  procedure  was  based  on  previous  work  by  Gurtin 
[71,731  Mikhlin  [98]  and  Sandhu  and  his  co-workers  [135-1381  The  self-adjoint  system 
of  equations  for  the  problem  is: 


A(u)  -  f  on  R  x [0 , oo ) 

Here  A  is  a  nutria;  of  operators.  Explicitly 


L(a)  p(,)/f+i«(i/k)  -t«-i-  o 


(i) 


-L 


**£ 

0 

0 

0 


0 

0 

-t* 


0 

-? 


0 

0 


-t* 

p 


0 

0 

-t* 

0 

t*oM8 


0  foMl.  t*M 


u 


where 


JL  t*(8  -3-  ♦  8 
2  '•  dh  »•  ai 


and 


p  - 

Also  in  (4A9X 


u. 

F 

m 

w 

■ 

G 

■ 

ir 

and  f  » 

0 

r-J 

0 

0 

i 

0 

(4^9) 


(4.70) 


(4.71) 


(4.72) 


(4.73) 
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Elements  of  the  operator  matrix  A  satisfy  self -ad  join  tneas  with  respect  to  the  bilinear 
mapping 


<f  ,g> 


f*gdR 


(4.74) 


La,  hey  satisfy  the  relation 

a  a 

Z>)-A,V*  -  <V  Za.jV«  *  D*(VUJ)  •  imUX . n 

r-i  r> 


(4.75) 


where  D^ti^.Uj)  denotes  quantities  asmciated  with  the  boundary  #R  of  the  region  of 
interest  R.  Consistent  boundary  conditions  for  (4j69)  are 
-t^u^j  «  -t^u-n^  00  S,x{0,oo) 


-t*  w[n1  ■  -t*WjOt 

on  SjX(0,oo) 

t*  trni  «  t*  inii  on 

Sj  x  (0 .  oo) 

(4.76) 

f  V,  -  fT.  on 

S4  x  (0 ,  oo) 

Consistent  form  of  the  internal  jump  discontinuities  is 

-^(u^j)'  «  -t'Cg,),^ 

on  S,(  x(0,oo) 

■ffwBj  ■  -t*gj 

on  Sa  x  [o ,  oo ) 

tMirn)’  «  t*gjnt 

on  S,(  x  (0 ,  oo) 

(4.77) 

<*( V/  *  t**4ni 

on  S4(  x  (0 ,  oo) 

Boundary  operators  C,,.  i. j  * 

1.2 . n  are  said  to  be  consistent  with 

the  matrix  of 

field  operators  if  in  (4.75) 


Dia(ui*uj)  *  <u**LcyV as  “ 


(4.78) 


r  i  n 

Here,  surfaces  S^,  S,,  S,  and  S,  are  embedded  in  the  interior  of  R  Operators  in  the 

self-adjoint  operator  matrix  in  (4.70)  have  the  following  relationships 

<t*u  ,r  =  —  <t*u  ,r  >_ 

'•j  'j  «  >*  <j.j  * 
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♦^iVVs,  + 

+  <f  (U.B/ .  r(j>Sj(  +  <fu, .  (rt/ n/>S4 


(4.79) 


•  #j>|  *  ""  <t*Wj  | » 


+  <t*wlH|  •  ^Sj  +  <*^,-wlS>*J 


(4.80) 


Here  we  unoe  that  <  .  >,  is  the  nun  of  quantities  evil  at  ted  over  the  subregions  of 
R  such  thst  all  the  surfaces  S,  S,,  $,  ut  contained  in  the  union  of  the  bound¬ 
aries  of  these  subregions.  For  the  coupled  system  (4^9)  the  governing  function  is 


defined  as; 


<i)  . 

0(u)  ■  <PU(.U>B  +  “  <*%.,.  V«  +  +  7^W.'W.>« 

+  ~  2<t*f.ir>t  +  <t*u, t).rtJ>K 

—  2<t*e(J ,  r(i>g  +  <t*(Eyu  +  or  MS^S^)  eu ,  e^>t 
+  2<t*otM8ifiJ,(>|  +  <t*Mf , £>t  -  2<ul,Fl>B  -  2<w. ,G.>b 
-<V*u|-2ul)nj>sl  "  <*.t*w|-2w,)n|>si 
+  <w(.t^w-2rr)ni>Sj  +  <u, .  t*(r p. -  2T,)  >s< 
-<rlJ,t^(uinJ)'-2(g1)1nJ)>Sti  -  <irIt^(wtn,)-2gJ)>S2i 
+  <w(.K(im|)'-2gJni)>Sji  +  <ui,f((riJnJ)'-2g4n1)>S4( 

The  Gateaux  differential  of  this  function  along  v  =  (u,,  w,,  W,  e,f  £}  is; 


(4.81) 


AyfKu)  =  <u,,pu.  +  p<a)w.-t*ry  ,j-2F.>b 

+  <u|,pu,  +  p<2>w,-t«T|JtJ>R 

(2)  . 

+  <W.  ,p(2)u  +(£—  +  1*— )w  -t*7rsib,i-2G  >„ 
||f  k  ‘  1  11 
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+  <  w,  .pl%  +  {dl  +  1*1)^,  -  fir  >, 

*  <Vt*ui.j”t%>a  *  <Vt*5».J_tV« 

+  <iy.-t^|J  +  KE,^  +  o2M«lj8u)eu  +  oM8^>t 

♦  <eij  •  **% +  +  a*M4IJV«u  +  *  MV>* 

+  <£,-t*ir  +  +t*aM8ueu  +  t^f>|1 

+  <£,-t*v  +  + 1*  or  +  t*M|  >  g 

-<r|JfKu^r2u^J)>s|  ~  <T,j»t*“inj>sl 

~<».t*(win1-2win1)>Si  -  <*  ,t*wtn(>Sj 

-<wi.f(tnxi-2^nl)>Sj  -  <w1,t^n,>Sj 

-<ul.^rlJnr2t,)>S4-  <Vt*Tijnj>s4 
-<?)i,t^(u1n/-2(g|)1nJ)>Sii  - 
-<w.tH(win1)'-2gP>Sji  -  < ir, fCwn,)' >Sa 
-<w1,f((imi)'-2gJnl)>  -  <wi,tK¥n,)'> 

^  —  »^...  .41  41 

-<3i,f((TlJn/-2Cg4)np>S4i  -  <u|,Krljn/>S4) 


(4.82) 


Using  (4.79)  and  (4.80)  the  Gateaux  differential  is  seen  to  vanish  if  and  only  if  all 
the  field  equations  along  with  the  boundary  conditions  and  the  jump  discontinuities  are 
satisfied.  (4.79)  and  (4.80)  relate  pairs  of  operators  in  the  operator  matrix  A  and  may 
be  used  to  eliminate  either  of  the  elements  in  each  pair  from  the  function  ft(u)  in 
(4.81).  Eight  alternate  forms  can  be  obtained  by  using  either  or  both  relations.  Elimina¬ 
tion  of  an  operator  A1J  from  the  function  implies  that  state  variable  Uj  need  not  be  in 
the  domain  M,}  of  A1J.  This  results  in  relaxating  the  requirement  of  smoothness  of  u^ 
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thereby  extending  the  spece  of  admissible  states.  In  the  context  of  the  finite  element 
method,  it  is  deer  that  the  extension  of  the  admissible  space  provides  greater  freedom 
in  selection  of  approximation  functions. 

If  the  admissible  state  is  constrained  to  satisfy  some  field  equation  and/or  bound¬ 
ary  conditions,  certain  forms  of  the  variational  principle  are  realized.  This 

procedure  is  used  to  reduce  the  number  of  free  variables  in  the  governing  function. 
Also,  certain  assumptions  in  the  spatial  or  temporal  variation  of  some  of  the  variables 
lead  to  approximate  theories.  In  the  context  of  direct  methods  of  approximation  the 
constarints  asm  mad  in  the  specialization  must  be  satisfied  by  admissible  states. 

As  an  example,  for  the  extended  functional  obtained  by  using  (4.79)  and  (4.80)  to 
eliminate  ru  ,  and  w ,  from  (4.81).  specialization  to  satisfy  (4.69),  and  (4.69)*  Le,  satis¬ 
fying  identically  the  kinematic  relationships  gives 

fil  J2>  1 

Q,  «  <pui,u1>J|  +  2 <p  Wj , Uj>j  +  <(^-  +  l*i-)w1.w1>J| 

+  <t*(Eliil  +  a2M81J5kJ)ekJ.e)J>B  +  2<t»or  >„ 

+  <t*M£.£>B  -  2<u,,Fj>B  -  2<wi,G,>B 
-2<T1J.t*(u1-at)nJ>Sj  -  2<ir.t^wi-wl)nl>S2 

-2<wi,t^n1>S3  -  2<u1,t*ti>s< 

—  2<r1J.t*((u1np'  — (g1)1np>Sjj  -  2<7T,t<(w.n1)'-g2)>S2i 

-  2< Wi * > s3|  -  2 <’ “l  ■  t’*4nl > s41  (4-83) 

If  the  field  variables  over  the  domain  are  continuous,  the  jump  discontinuity  terms 
drop  out  giving  the  specialization; 

ft,  =  <Pui.u,>B  +  2<p(2)w.,u.>R  +  <(£r  +  l»l)w..wi>B 
+  <t*(Ejjk|  +  a  ,eij>R  +  2<t*aM81.ejj,f>| 
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+  -  2<ui.F1>I  -  2<w(.Gl>E 

-2<fjJ1Ktti-u|)nJ>Ji  -  2<*,Kwl-w()n|>Si 

—  2<w1,t*ih»,>t^  —  2<u(,t*fj>g^  (444) 

Further  of  (444)  to  the  case  when  displacement  boundary  conditions  an 

identically  satisfied  yields  the  function  governing  the  two  field  formulation  proposed 
by  Ghabouan  [S3]  except  that  in  the  present  formulation  the  boundary  terms  an  con¬ 
sistent 

Figure  7  diagrams tically  depicts  the  possible  extensions  of  the  general  variational 
principle  based  on  the  direct  formulation.  Figure  8  shows  the  same  for  the  comple¬ 
mentary  formulation.  In  either  case,  only  the  specialisations  listed  in  the  report  (item 
1.7  in  Appendix  B)  an  shown.  Evidently,  other  extended  forms  could  be  used  as  start¬ 
ing  points  for  specialization. 

Details  of  this  effort  an  contained  in  items  1.7  and  34  in  Appendix  B.  For  non¬ 
linear  problems,  a  quaailin  earned  form  of  the  nonlinear  equations  was  used  to  develop 
a  variational  formulation  (item  1.13  in  Appendix  B). 

4.4.2  Two-  and  Three-Field  Formulations 

In  the  two-field  formulations  of  Biot's  theory  the  soil  displacement  and  the  dis¬ 
placement  of  fluid  relative  to  the  soil  wen  used  as  the  two  field  variables.  The  pore 
pressures  wen  determined  through  a  constitutive  relationship  using  volumetric  strain  of 
the  soil  and  the  change  in  fluid  content.  During  the  course  of  the  present  research,  it 
was  felt  that  this  approach  may  not  yield  sufficiently  accurate  estimates  of  the  pore 
pressures  because  of  the  need  for  numerical  evaluation  of  the  derivatives.  A  three-field 
formulation  introducing  the  pore-pressure  as  the  third  field  variable  was  derived  as  a 
specialization  from  the  general  variational  principle.  Asstuning  the  boundary  conditions 


n(u) 

(Governing  Functional) 


Fi£urt  7:  Family  of  Variational  Principles  by  Direct  Formulation 
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on  Sj  and  S,  are  identically  satisfied,  the  governing  function  for  this  was: 


a  -  <pui,ul>8  +  2<p(2)wj,ui>8  +  <(p(2)/f  + 

-2<t*iri,wi>R  +  <t*Eijkieu,e.j>R  +  2<t*a«ijelj,ir>R 
-  <t*Z.,ir>R  -  2<ui,Fi>R  -  2<wl.Gi>R 

+  2<ir,t*w.n,>c  —  2<u  ,t*f>.  (4.85) 

1  1  *  5^ 

This  formulation  was  implemented  in  a  finite  element  computer  code  to  obtain  continu¬ 
ous  pore  pressure  distributions.  The  three-field  formulation  also  allowed  direct  specifica¬ 
tion  of  the  fluid  pressures  on  the  boundaries,  which  is  not  possible  with  the  two-field 
formulations,  where  this  boundary  condition  could  only  be  applied  as  a  linear  con¬ 
straint  in  terms  of  soil  displacement  and  relative  fluid  displacement.  The  studies 
showed  that  though  the  numerical  difference  in  the  results  from  the  two-  and  three- 
field  formulations  was  only  slight,  the  three-field  formulation  was  much  more  expen¬ 
sive.  However,  it  has  the  distinct  advantage  of  being  able  to  prescribe  boundary  values 
for  pore  pressures.  Items  1.10,  3.4  and  4.1  of  Appendix  B  contain  details  of  the  for¬ 
mulation  and  a  study  of  its  effectiveness. 

4.4.3  Spatial  Discretization 

Most  succesful  schemes  for  approximate  solution  of  the  coupled  problem  of  quasi¬ 
static  soil  deformation  and  fluid  flow  have  been  based  on  the  use  of  higher  order 
interpolation  for  the  displacements  and  a  lower  order  for  the  pore  water  pressures. 
These  elements  are  expensive  to  use  in  terms  of  computer  time.  Elements  based  on  the 
use  of  the  same  order  of  interpolation  had  been  found  to  give  unreliable  results  just 
after  loading  and  would  be  suspect  for  use  in  dynamic  problems.  Incompatible  interele¬ 
ment  boundaries  have  been  used  to  combine  the  economy  of  lower  order  elements  and 
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yet  retain  higher  order  local  interpolation.  These  elements  satisfy  the  patch  test  for 
completeness  for  certain  element  geometries.  Ghaboussi's  four  point  incompatible  element 
[52^3]  was  implemented  in  a  computer  program  for  consolidation  analysis  to  compare 
its  performance  with  that  of  higher  order  elements  for  efficiency  and  accuracy.  Details 
of  this  study  are  contained  in  item  13  of  appendix  B.  These  studies  were  useful  in 
the  selection  of  the  strategy  appropriate  to  the  problem  of  dynamic  analysis  of  liquid- 
filled  soils.  The  comparative  study  showed  that  the  incompatible  element  gave  results 
almost  identical  to  those  obtained  using  the  higher  order  elements  based  on  biquadratic 
interpolation  for  displacements,  but  was  significantly  more  economical.  This  made  Gha¬ 
boussi's  element  a  good  candidate  for  extension  to  nonlinear,  three-dimensional  and 
dynamic  problems. 

Three  different  finite  element  strategies  were  used  to  cover  the  cases  of  one-  and 
two-dimensional  wave  propagation.  Both  bilinear  and  biquadratic  interpolation  schemes 
were  implemented  along  with  Ghaboussi's  incompatible  element  and  cubic  Hermite  poly¬ 
nomials  for  one-dimensional  wave  propagation.  These  interpolation  schemes  were  imple¬ 
mented  in  a  dynamic  analysis  computer  program  and  used  to  solve  one-dimensional 
wave  propagation  problems  for  which  exact  solutions  were  available.  The  code  was 
used  to  solve  one-dimensional  wave  propagation  in  a  single  material.  Both  the  steady 
state  and  the  transient  cases  were  considered.  Garg's  [50]  theoretical  solutions  for  the 
case  of  weak  as  well  as  strong  couplings  were  used  as  benchmarks.  Application  of 
the  computer  code  to  all  these  problems  showed  excellent  agreement  between  the 
numerical  and  the  theoretical  solutions.  The  code  could  not  be  tested  for  two- 
dimensional  wave  propagation  because  of  lack  of  exact  solutions  for  that  case.  Reports 
listed  as  items  1.3,  1.15  and  23  of  Appendix  B  describe  some  of  this  work. 
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AAA  Singularity  Elements 

Conventional  finite  element  procedures  cannot  adequately  model  pore  pressures 
near  loaded  free-draining  boundaries  because  of  the  singularity  in  the  pore  pressure  that 
exists  immediately  after  the  load  is  applied.  Special  finite  elements,  developed  by  Lee 
[89l  to  allow  for  line  singularities  were  used  to  simulate  propagation  of  waves  in  sat¬ 
urated  soil.  These  elements  use  interpolating  functions  of  the  type 

1  —  ax  —  (1  —  a)x° 

where  the  index  n  is  sufficiently  large  and  coefficient  a  is  chosen  to  be 
a  =  1  —  exp(—  mt) 

and  were  found  to  be  satisfactory  for  proper  representation  of  singularity  at  the  wave 
front.  These  elements  have  the  property  that  for  t  =  0,  a  =  0  and,  therefore,  the  interpo¬ 
lating  function  is  1— x“  and  as  t  increases  the  function  approaches  1  —  i.  This  element 
reproduces  the  line  singularity  occuring  in  one-dimensional  consolidation  problem  imme¬ 
diately  after  loading  and  the  singularity  at  a  wave  front.  The  interpolating  scheme 
would  approach  1  —  x  as  time  t  increases.  This  research  is  described  in  items  1.1  and 
3.1  in  Appendix  B. 

AAJS  Time  Domain  Integration 

The  results  for  one-dimensional  analysis,  for  which  exact  solutions  are  available, 
showed  that  the  conventional  time-domain  integration  procedures  found  to  be  quite 
effective  for  single  material  problems,  were  also  acceptable  for  the  coupled  problem  of 
dynamics  of  saturated  soils.  A  popular  scheme  is  Wilson's  0,  y,  0  single  step  method. 
For  two-dimensional  cases,  no  exact  solutions  are  available.  For  this  reason,  all  the  code 
verification  had  to  be  done  on  one-dimensional  wave  propagation  problems.  The 
requirements  for  an  acceptable  time  step  integration  scheme  were  that  the  results 
should  be  insensitive  to  the  choice  of  time-domain  integration  selected  by  the  user  and 
that  with  reduction  of  the  size  of  the  time  step  the  approximate  solution  should  con- 
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verge. 


A  modification  to  Wilson's  method  was  introduced  in  order  to  specify  velocity 


boundary  conditions.  Finite  element  discretization  of  the  spatial  domain  leads  to: 
MU  +  CU  +  KU  =  R 


(4.86) 


where  M,  K  and  C  axe  the  mass,  stiffness  and  damping  matrices,  respectively,  and  R  is 
the  load  vector  at  time  tB  In  (4.86),  the  square  brackets  and  curly  braces  have  been 
dropped  for  convenience,  as  in  (4.52).  In  Wilson's  scheme  the  displacements  and  veloci¬ 
ties  at  time  (t  +  0At)  are  expressed  in  terms  of  U,  U  and  U  at  time  t„  as 


Un+a  *  u„  +  #AtUn  +  (1/2  —  0)(0At)2U  +  0C0At)2l) 


(4.87) 


Un+0  =  U  +  (l-y)(0At)U 


(4.88) 


in  which  0  and  y  are  Newmark's  coefficients.  Substituting  (4.87)  and  (4.88)  into 
(4.86)  yields 


KU  .  =  R  - 

mw  irv 


(4.89) 


where 


K  =  K  +  - - — -M  +  SC 

0(0At)2  W 


Rf^  +  +  Wm)03^ 


(4.90) 


a^  =  Uo  +  (0At)Un  +  (l/2-0)(0At)2Ut 
=  U  +(l-0/y)(0At)2Un 


Assuming  cubic  variation  of  nodal  point  displacement  over  the  time  step  (tB ,  %  + ,)  in 
terms  of  displacement,  velocity  and  acceleration  at  time  tB,  the  values  of  these  quanti¬ 


ties  at  time  (tB  +  At)  are 


u~.  -  >-4u.+  '-ir^.  +  K'-O'402'3. 
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L.  =  y  Tcu^-un)  + 

003At2 


1  - 


002 


k  +  0 


u  =  1 Cu.-u  ) - - — u  +(i — —  \j 

-*1  0038t2  “  002 At  ■  V  ° 


200 
20? 


(4.91) 


In  a  given  problem,  when  velocities  are  prescribed  at  certain  points,  the  corre¬ 
sponding  accelerations  and  displacements  are  also  known.  In  developing  the  modified 
scheme,  all  features  of  Wilson's  method  were  sought  to  be  retained.  This  enabled  appli¬ 
cation  of  the  scheme  directly  without  elimination  of  known  degrees-of-freedom  which 
is  quite  cumbersome  in  dynamic  problems.  Let  subscripts  a,  b  denote  the  unknown  and 
specified  quantities,  respectively.  Then  for  the  stage  (n  +  0)  at  time  tn  +  0At,  (4.86)  can 
be  rewritten  in  partitioned  form  as 


M,b  |u 

K.  fc| 


(a*4) 


-  c 

"«»  '-'»b 

-  c 

■'b*  '■'bb 


(a*4) 


CJK 

h>b  “b 


(a*9) 


Rb 


(n+fl) 


Rearrangement  of  terms  in  (4.92)  gives 


(4.92) 


(4.93) 


From  (4.87) 

=  0(0AO2'  [U“^  ~  U°  ”  9 At  ““  ~il/2~  ® (0 At)2“  J 
Substituting  (4.94)  in  (4.88) 

Rewriting  (4.94)  and  (4.95)  for  the  unknown  quantities  u,  and  ut 


J  _  1 

U. 

_ 

h- 

-  (0At) 

K 

-  (1/2  —  0)(0At)2 

0(0At)2 

(a*«) 

0 

(n*4) 

k> 

a 

0 

D 

0 

Q 

u 

u 

u 

,(u 

0 

0(0At) 

B 

— 

o* 

-  (0At)(l  —  0/y) 

0 

—  (1/2  — 0/y)(0At)2r 

(n+0) 

n 

n 

n 

(4.94) 


(4.95) 


(4.96) 


(4.97) 


Ill 


i 


Substituting  (4.96)  and  (4.97)  into  (4.93)  gives 


kIJki 


(4.98) 


where 


K“  =  +  +  W£TC‘ 


K" =  K“  +  + 

+  +  ^foCbb 


(4.99) 


(4.100) 


(4.101) 


(4.102) 


:  =  *  [m] — - — 

!»  W,  -  ««4t) 

0  l,a*6) 


1  U  (O  U 

—2 ■  *  +  „  +  (9At)  * 

At)2  (0  M  |0) 

a  o*6  i 


+  WAon  +  (1/2  —  0)  (0At)2|U*' 


+lclsfefcj +  P 

n  i 

-h°j  -ny 

-  ■  A  '  '  -  .fl 


+  (0At)(l-0/y*  ■)  +  (l/2-/3/y)(0At)2 


(4.103) 


Wilson's  method  of  allowing  for  prescribed  displacement  boundary  condition  was  used 
to  rewrite  (4.98)  as 


K»  0 K 

O  I  “b 


_  R. 


(4.104) 


where 


R*  =  R«  K*b^W1^) 


(4.105) 
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and  ub  are  the  prescribed  values  of  ub  I  denotes  the  identity  matrix.  The  first  9et  of 
equations  in  (4.104)  is  the  first  set  of  (4.98)  modified  for  the  known  ub  and  the  sec¬ 
ond  set  of  equations  are  trivial  equations  introduced  to  avoid  reordering  of  elements  of 
the  displacement  vector. 


The  procedure  discussed  above  was  used  to  solve  the  problem  of  wave  propaga¬ 
tion  in  a  finite  soil  column  of  length  (L  =  50cm.)  shown  in  Figure  9.  The  excitation 
was  applied  at  the  top  (x  -  0).  The  base  of  the  soil  column  was  assumed  to  be  rigid 
and  impervious.  The  material  properties  were  chosen  to  be  the  same  as  used  by  Garg 
[50l  Le, 

n(1>  =  0.82,  n<2)  =  0.18,  p(1)  =  2.1812  g/cm3,  p(2)  =  0.18  g/cm3 
Kj  =0.36  x  1012,  K  =  0.118  X  1012,  K2  =0.22  x  1011,  G  =  0.99xl0u  (all  in  dynes/cm2) 
These  are  related  to  Biot's  constants  and  yield 

E  =  0.2321  x  1012  dynes/cm2,  v  =  0.171,  a  =  0.6722,  M  =  0.1047  X  1012  dynes/cm2 
Two  example  problems  with  different  boundary  conditions.  Figure  10,  one  following 
Garg  [50]  and  the  other  suggested  by  Morland  (item  1.15.  in  Appendix  B)  were  solved 
using  the  numerical  inversion  of  the  Laplace  transform  solution  and  the  direct  finite 
element  procedure.  Mathematically,  the  boundary  conditions  are: 


Example  1: 

v(1)(0,t)  =  H(t) 


(4.106) 


v<2)(0 ,  t)  = 
where  H(t)  is  the 


H(t)  j 

Heaviside  function.  Corresponding  velocity  transformations  are 


v^O.p)  =  1/p 
v2)(0,p)  =  1/p 


(4.107) 
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9:  Representative  Soil  Column 
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Figure  10:  Types  of  Excitation  Applied  on  the  Soil  Column 
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Example  2: 


v(1)(0,t)  =  H(t) 


(4.108) 


v(2)(0,t)  =  [1  — 0.2e~t/T]H(t)  ] 

t  =  (L/C,)  =  140.9  fistc  was  taken  as  a  normalizing  factor.  In  this  example,  the  fluid 
velocity  specified  at  the  boundary  was  different  from  the  specified  solid  velocity  and 
increased  gradually  from  a  value  of  0.8  at  t  =  0  to  unity  over  the  time  scale.  The 
velocity  transformations  are  given  by 

^°(0,p)  =  1/p  1 


(4.109) 


v^O.p)  =  1/p  — 0.2/(p  +  l)  | 

Two  values  of  the  drag  parameter 

D  =  0.219  x  102  g/cm3— sec 

D  =  0.219  x  10*  g/cm3-sec 

representing  the  so-called  low  drag  (free  relative  motion  between  constituents)  and  high 
drag  (negligible  relative  motion)  were  used.  The  corresponding  values  for  the  ratio 
K/n  used  in  the  finite  element  analysis  were  0.148  x  10  J  and  0.148  x  10'*,  respectively. 


Numerical  inversion  of  the  Laplace  transform  solution  to  the  problem  was  carried 
out  using  5000  terms  in  Dubner's  [41]  formula 

oo 

f(t)  =  (l/fOe”  Re{f(r  +  k7ri/2r)  cos  (kirt/2T)]  0<t<r  (4.110) 

k-0 

where  prime  signifies  that  only  half  of  k  =  0  term  is  included  in  the  sum.  Dubner 
[41]  showed  that  the  error  could  be  made  small  for  0<t<r  by  choosing  rr  sufficient¬ 
ly  large,  where  r  is  a  real  number.  The  value  of  rr  was  set  equal  to  5.017.  The 
velocity  histories  at  four  locations,  namely  10  cm  and  30  cm  from  the  free  surface. 
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were  recorded  over  a  total  time  period  of  986.4  /usee  at  intervals  of  14  ptsec.  This 
allowed  for  six  reflections  of  the  faster  C+  wave  and  two  reflections  of  the  slower  C_ 

wave.  The  CPU  time  on  the  IBM  3081  mainframe  was  70  min. 

The  finite  element  solution  was  a  uniform  spatial  mesh  for  Example  1  but  for 
Example  2  modelled  the  boundary  layer  by  a  fine  mesh  of  elements  near  the  top  sur¬ 
face  including  a  singularity  element  adjacent  to  the  surface.  Thus  in  Example  1  the 
spatial  discretization  consisted  of  100  linear  elements  and  986  time  steps  of  size  1  /usee 
in  the  time  domain.  In  Example  2,  100  elements  Of  0.005  cm  length  were  employed 
near  the  top  surface  and  100  elements  of  0495  cm  for  the  remainder  of  the  column. 
The  temporal  integration  involved  141  steps  of  0.01  /usee  and  985  steps  of  1  //sec. 

The  velocity  histories  for  low  drag  and  high  drag  for  Example  1  are  shown  in 
Figure  11  to  Figure  14  and  Figure  15  to  Figure  18,  respectively.  For  Example  2,  Figure 
19  to  Figure  22  illustrate  the  low  drag  effects  and  Figure  23  to  Figure  26  the  high 
drag.  In  both  examples,  excellent  agreement  of  numerical  results  is  seen.  In  example 
2,  a  refined  mesh  consisting  only  of  linear  elements  was  not  able  to  reproduce  the 

sharp  wave  fronts  and  large  oscillations  were  encountered.  This  was  overcome  by  the 

use  of  singularity  element  near  the  top  boundary.  The  shape  functions  for  this  element 
were  of  the  type  (!—{*)  and  over  (0,1).  The  index  0  was  taken  to  be  100. 
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Figure  13:  Example  1:  Solid  Velocity  History  at  30  cm.  (Low  Drag) 
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Figure  14:  Example  1:  Fluid  Velocity  History  at  30  cm.  (Low  Drag) 
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Figure  16:  Example  1:  Fluid  Velocity  History  at  10  cm.  (High  Drag) 
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f  ij>ure  l7:  Example  1:  Solid  Velocity  History  at  30  cm.  (High  Drag) 
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Figure  J8:  Example  1:  Fluid  Velocity  History  at  30  cm.  (High  Dtag) 
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Figure  20:  Example  2:  Fluid  Velocity  History  at  10  cm.  (Low  Drag) 


Figure  22:  Example  2:  Fluid 


f  igure  24:  Example  1 
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4X6  Nonlinear  Problems 

For  nonlinear  problems,  an  incremental  approach  was  necessary.  The  equations 
governing  this  case  along  with  a  variational  formulation  of  the  problem  were  devel¬ 
oped  (item  1.13  of  Appendix  B).  Only  material  nonlinearity  was  considered. 

At  any  instant  of  time  t,  the  equilibrium  forces  acting  on  the  discretized  system 
can  be  represented  by 


Mh  M„ 
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0  0 

u  K  Ou 

1  (u 

F| 

P. 

•  • 

4- 

oc„ 
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Alternatively,  one  might  write 

F,(t)  4-  FD(t)  +  F*(t)  +  F“(t)  =  P(t)  (4.112) 

where  Fp  F0  are,  respectively,  the  inertial  and  the  damping  force  vector.  Fj  represents 
the  internal  resisting  force  related  to  the  9olid  deformation  only  and  F^  is  the  internal 
resisting  force  arising  out  of  deformation  of  the  fluid  and  coupling  between  the  two 
phases,  P  denotes  the  applied  load  vector.  A  short  time  At  later,  the  equation  would 
be 

F,(t  4-  At)  4-  FD(t  +  At)  +  F*(t  +  At)  4-  F~(t  +  At)  =  P(t  4-  At)  (4.113) 

Subtracting  (4.112)  from  (4.113), 

F,(t  +  At)  -  F,(t)  +  FD(t  4-  At)  -  FD(t)  4-  F*(t  4-  At)  -  F*(t) 

4-  F~(t  4-  At)  -  F*(t)  =  P(t  4-  At)  -  Kt)  (4.114) 

Noting 

F,(t  4-  At)  *  M(t  4-  At)u(t  4-  At),  (4.115) 

expanding  M(t  4-  At)  and  u(t  4-  At)  in  terms  of  Taylor's  series,  and  retaining  only  the 
first  order  terms  in  At  gives 

Fj(t  4-  At)  =  M(t)  u(t)  4-  M(t)Au(t)  4-  AM(t)\i(t)  (4.116) 


Similarly,  the  quantities 
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FD(t  +  At)  =  C(t  + At)u(t) 

F*(t  +  At)  =  K*(t  + At)u(t  + At)  (4.117) 

F*(t  +  At)  =  K**(t  +  At)  u(t  +  At) 

can  be  approximated  as 

FD(t  +  At)  =  CXt)u(t)  +  CXt)Au(t)  +  A(Xt)u(t) 

F*(t  +  At)  =  K*(t)u(t)  +  K*(t)  Au(t)  +  AK*(t)  u(t)  (4.118) 

F~(t  +  At)  =  K"(t)u(t)  +  K"(t)  Au(t)  +  AK*  u(t) 

Use  of  (4.116)  and  (4.118)  in  (4.114)  yields 

M(t)Au(t)  +  AM(t)u(t)  +  C(t)Au(t)  +  AC(t)u(t)  +  K*(t)Au(t)  +  AK*(t)u(t) 

+  K“  Au(t)  +  AK*  u(t)  =  Kt  + At)  -  M(t)u(t)  +  CXt)u(t) 

—  K*(t)u(t)  +  K“(t)u(t)  (4.119) 

This  represents  a  general  form  of  incremental  equations.  If  mass,  damping  and  stiffness 
quantities  at  time  t  are  known,  (4.119)  can  be  solved  for  Au(t)  by  step-forward  inte¬ 
gration  scheme,  which  also  yields  Au(t)  and  Au(t).  In  doing  so,  the  quantities  them¬ 
selves  are  dependent  on  the  solution  Au(t)  and  hence  an  iterative  scheme  to  reduce  the 
cumulative  error  is  necessary. 

The  theory  was  specialized  to  the  case  of  nonlinearity  only  in  the  soil/stress  rela¬ 
tions.  This  case  was  implemented  in  two  finite  element  programs.  In  the  code  NAOWP 
(Nonlinear  Analysis  of  Wave  Propagation)  elastic- perfectly  plastic  and  bilinear  stress- 
strain  relations  were  used.  The  other  code  named  DANS  (Dynamic  Analysis  of  Nonli¬ 
near  Soils)  incorporated  a  more  general  model  of  elastic-plastic  work-hardening  proposed 
by  Singh  [160].  A  modular  structure  was  used  so  that  a  variety  of  models  could  be 
selected.  Local  iteration  was  employed  within  each  time  step  and  convergence  assured 
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before  going  on  to  the  next  step.  NAOWP  was  designed  to  model  one-dimensional  wave 
propagation.  linear  variation  over  the  spatial  elements  was  assumed.  In  DANS  (item 
1.19  of  Appendix  B),  bilinear  and  biquadratic  isoparametric  elements  were  implemented 
for  two-dimensional  wave  propagation  analyses.  Since  no  exact  solutions  for  nonlinear 
wave  propagation  in  fluid-saturated  soils  were  available,  the  codes  were  verified  against 
exact  solutions  for  single  material  wave  propagation.  Exact  solutions  for  bilinear  solids 
subjected  to  dynamic  excitation  have  been  developed  by  Belytschko  [10]  and  for  an 
elastic-perfectly  plastic  solid  by  Wood  [173].  Figure  27  describes  the  discretization  for 
a  soil  column  as  well  as  the  suddenly  applied  loading  for  the  three  cases  tested.  Figure 
28  shows  the  stress  history  for  case  3  plotted  against  the  exact  solution  by  Belytschko. 
Figure  29  shows  the  stress  pulse  of  short  duration  used  by  Wood  [173]  along  with  the 
elastic-perfectly  plastic  soil  column.  Figure  30  to  Figure  34  show  the  stress  profiles  at 
time  equal  to  4,  8,  12,  16  and  20  plotted  against  the  exact  solution  by  Wood  [173]  for 
different  levels  of  mesh  refinement.  Item  1.13  of  Appendix  B  contains  details  of  the 
approach  as  well  as  illustrative  applications  of  the  two  computer  programs  to  wave 
propagation  through  a  saturated  soil  layer. 


Laterally  Restrained 


29:  Solid  Column  Under  Stress  Pulse  of  Short  Duration 
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Figure  34:  Stress  Profile  at  t 


Section  V 


LABORATORY  INVESTIGATIONS 


5.1  INTRODUCTION 

Historically,  laboratory  studies  of  material  behavior  of  saturated  sands  under  cycl¬ 
ic  stress  conditions  have  played  an  important  role  in  the  determination  of  liquefaction 
potential.  Numerous  investigators  have  tried  to  model  and  predict  the  potential  and 
probability  of  liquefaction  occurring  in  soils.  Various  test  apparatus  have  been  designed 
or  modified  in  an  attempt  to  provide  an  accurate  representation  of  the  stress  state  gen¬ 
erated  in-situ  by  ground  motion.  A  number  of  experimental  devices  including  the  cyclic 
triaxial,  cyclic  simple  shear,  torsional  shear  and  shaking  table  have  been  developed.  A 
wide  diversity  of  data  have  been  generated  with  the  use  of  a  dynamic  (or  pseudo¬ 
dynamic)  excitation  loading  pattern.  Detailed  reviews  of  these  experimental  programs 
and  design  methods  based  upon  them,  can  be  found  in  several  state-of-the-art  reports. 
In  particular,  those  by  Seed  [155-1581  Finn  [44l  Casagrande  [31]  and  the  National 
Research  Council  [109]  should  be  noted. 

The  laboratory  method  used  in  this  study  to  load  the  sand  samples  was  a  shak¬ 
ing  table.  It  is  well  known  that  the  variety  of  small  scale  apparatus  currently  in  use 
in  the  laboratory  introduce  non-uniform  stress  and  strain  fields  in  the  sample  being 
tested  (eg.  [31,127,94,174]  ).  Since  the  onset  of  liquefaction  is  clearly  a  local  phenom¬ 
enon,  it  is  certain  that  a  measure  of  the  liquefaction  resistance  of  a  saturated  sand 
would  be  affected  by  these  stress  and  strain  concentrations.  The  result  of  stress  concen¬ 
tration  induced  liquefaction  would  be  an  underprediction  of  the  true  liquefaction  resis- 


tance.  Therefore,  although  the  shaking  table  is  not  as  widely  used  as  the  other  labora¬ 
tory  devices  identified  above,  because  it  does  more  closely  simulate  actual  field  condi¬ 
tions,  several  previous  investigators  have  used  it  as  the  basis  for  liquefaction  studies. 
Among  the  many  reports  available  on  experimental  programs  are  several  in  which 
shaking  table  tests  were  conducted  specifically  for  the  purpose  of  studying  liquefaction 
including  the  programs  of  Finn  [4Sl  O-Hara  [llSl  DeAlba  [38],  Seed  [138]  and  Sasaki 
[139]. 


In  an  experimental  program  performed  on  samples  of  saturated  sand  in  1971, 
Finn  [43]  demonstrated  the  usefulness  of  the  shaking  table  for  conducting  liquefaction 
studies.  They  observed  that  the  shaking  table  offered  several  advantages  over  cyclic  tri- 
axial  and  simple  shear  devices.  Chief  among  these  advantages  were:  embedded  instru¬ 
mentation  having  a  negligible  effect  on  sample  response  can  be  used;  the  distribution  of 
pore  pressures  over  time  can  be  monitored;  the  uniform  accelerations  developed  in  the 
plane  strain  specimens  more  closely  corresponds  to  actual  field  conditions. 

O-Hara  [l  15]  conducted  shaking  table  tests  on  two  different  uniform  sands.  In 
addition,  he  performed  cyclic  triaxial  and  cyclic  simple  shear  tests  on  the  same  materi¬ 
al.  He  observed  that  sand  samples  tested  on  his  shaking  table  typically  showed  an 
increased  resistance  to  initial  liquefaction  when  compared  with  the  behavior  of  the 
same  materials  when  tested  in  either  of  the  two  small  scale  devices. 

The  shaking  table  tests  reported  by  DeAlba  [38]  performed  on  specimens  of  Mon¬ 
terey  No.  0  sand.  The  samples  were  4  inches  high  by  90  by  42  inches  at  the  base 
tapering  to  74  by  30  inches  at  the  top.  This  shape  was  chosen  so  that  a  rubber  mem¬ 
brane  could  easily  be  placed  over  the  specimen  and  then  pressurized  to  simulate  a  bur¬ 
ied  soil  element  The  size  was  chosen  to  provide  free  field  conditions  in  the  central 
portion  of  the  specimen.  The  membrane  under  which  the  sample  was  confined  allowed 
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the  specimen  to  deform  during  the  application  of  the  cyclic  load.  Additional  surcharge 
was  applied  to  the  specimen  in  the  form  of  a  reaction  mass  which  was  composed  of 
steel  shot  placed  in  a  bag  on  the  top  surface  of  the  sand  sample.  Pore  pressures  were 
measured  at  several  locations. 

The  results  of  further  tests  conducted  on  the  Berkeley  apparatus  were  presented 
by  Seed  [158]  in  an  experimental  program  designed  to  study  the  effects  of  seismic  his¬ 
tory  on  a  sand  deposit.  In  these  tests.  Seed  etal  observed  a  dramatic  increase  in  the 
number  of  cycles  required  to  induce  liquefaction  in  samples  which  previously  had  been 
subjected  to  cyclic  motions  significant  enough  to  raise  the  pore  pressure  but  not  enough 
by  themselves  to  cause  liquefaction.  Seed  etaL  attribute  the  observed  increase  in  cyclic 
strength  to  grain  rearrangement.  They  pointed  out  that  there  is  substantial  evidence 
that  these  higher  values  of  liquefaction  resistance  are  more  representative  of  the  actual 
performance  of  natural  sand  deposits  which  have  been  subjected  to  past  cyclic  motions. 

Sasaki  f  139]  conducted  a  series  of  shaking  table  tests  on  sand  samples.  The  sample 
size  they  chose  was  much  larger  than  had  been  used  in  any  previous  programs  (12 
meters  long  by  3  meters  high  by  2  meters  wide)  and  it  is  therefore  more  likely  that 
true  free  field  conditions  existed  in  their  samples.  The  samples  were  formed  by  pluvia- 
tion  through  air  similar  to  the  procedure  employed  by  Seed  [158].  Saturation  was 
reportedly  accomplished  by  displacing  the  air  through  the  infiltration  of  water  from 
the  bottom  of  the  sample.  The  sand  surface  was  unconfined.  Instrumentation  consisted 
of  embedded  pressure  transducers  and  accelerometers.  Measurements  of  cyclically 
induced  increases  in  pore  pressure  were  made  at  a  total  of  35  locations  within  the 
sample  in  both  the  free  field  and  in  the  vicinity  of  an  embedded  concrete  box  intend¬ 
ed  to  simulate  a  roadway.  It  should  be  noted  that  the  increased  pore  pressure  was  an 
order  of  magnitude  greater  immediately  beneath  the  roadway  than  it  was  in  the  free 


field.  Clearly  the  rate  of  pore  pressure  increase  during  shaking  is  very  sensitive  to 
the  presence  of  local  irregularities,  in  this  case  rigid  inclusions. 

In  order  to  minimize  these  local  effects  on  an  experimentally  obtained  estimate  of 
liquefaction  potential,  the  test  apparatus  described  in  the  following  sections  was 
employed.  Herein  we  give  a  summary.  Details  are  given  in  items  1.9,  2.6  and  4.3  of 
Appendix  B. 

5.2  EXPERIMENTAL  FACILITIES 

52.1  Introduction 

The  facilities  used  in  the  laboratory  investigation  of  the  liquefaction  phenomenon 
were  designed  for  the  purpose  of  providing  reproducible  results  in  which  the  stress  and 
strain  conditions  at  the  sample  boundaries  were  well  understood  and  could  be  recon¬ 
structed  accurately  in  a  numerical  model.  They  consisted  of  a  unidirectional  shaking 
table  to  which  was  attached  a  test  box  with  a  capacity  if  approximately  one  cu.  ft.  of 
soil.  The  instrumentation  employed  to  monitor  the  sample  behavior  consisted  of  accel¬ 
erometers  and  pressure  transducers. 

5.2~2  Shaking  table 

A  shaking  table  with  a  capacity  of  2500  lbs.  located  in  the  Soil  Dynamics  labo¬ 
ratory  at  the  Ohio  State  University  was  designed  for  the  liquefaction  testing  program. 
Table  motion  was  provided  by  an  MTS  system  capable  of  producing  peak  accelerations 
of  approximately  ±  2  g.  These  high  acceleration  levels  were  achieved  when  the  origi¬ 
nal  3  gpm  pump  was  replaced  by  a  10  gpm  pump  purchased  with  contract  funds.  In 
order  to  make  best  use  of  the  increased  capacity  of  the  pump,  the  existing  single  10 
gpm  servovalve  was  replaced  by  dual  10  gpm  servovalves  in  conjunction  with  two 
accumulators.  The  accumulators  significantly  improved  the  system  performance  during 


conditions  of  peak  fluid  flow.  The  input  signal  could  be  either  periodic  or  an  exter¬ 
nally  programmed  random  acceleration  time  history.  A  schematic  diagram  of  the  lique¬ 
faction  testing  system  is  shown  in  Figure  35. 

5.23  Test  Box 

The  test  chamber  which  was  bolted  to  the  shaking  table  was  designed  so  that: 

•  The  length  to  height,  ratio  of  the  samples  would  be  such  that  a  free  field  plane 
strain  condition  existed  in  a  substantial  portion  of  the  specimen. 

•  Moveable  inner  walls  which  would  provide  the  required  lateral  support  to  the 
sample  during  construction,  could  be  withdrawn  at  the  start  of  the  cyclic  test  so 
t^at  the  sample  could  deform  under  plane  strain  conditions. 

•  Normal  stresses  on  the  horizontal  faces  of  the  sample  could  be  chosen  indepen¬ 
dently  of  the  stresses  on  the  top  face,  thereby  allowing  for  tests  to  be  construct¬ 
ed  on  samples  without  requiring  that  K  be  -  1. 

A  diagram  of  the  test  chamber  is  shown  in  Figure  36. 
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Figure  35:  Large  Scale  Liquefaction  Testing  System 


ADJUSTING  BOLTS  FOR  INNER  WALLS 


5J-4  Rubber  Membrane 


The  rubber  membrane  surrounding  the  sample  was  made  to  the  exact  dimensions 
of  the  sample  test  chamber  with  the  inner  walls  in  their  unretracted  position.  It  was 
constructed  so  that  different  horizontal  and  vertical  pressures  could  be  applied  to  the 
sample.  The  membrane  had  to  be  thin  and  flexible  so  that  the  confining  pressures 
would  be  transmitted  uniformly,  but  strong  enough  to  withstand  the  working  pressures 
used.  A  fiber  reinforced  (nylon  on  nitrile)  synthetic  rubber  sheet  with  a  thickness  of 
.025  cm  met  these  requirements  and  was  to  make  the  membranes.  The  lower  mem¬ 
brane,  to  which  the  horizontal  pressure  was  applied,  was  in  the  shape  of  a  rectangular 
box,  supporting  the  bottom  and  four  sides  of  the  sample.  The  top  membrane  was  a  rec¬ 
tangular  sheet  placed  on  top  of  the  sample  after  construction.  With  the  lid  of  the  test 
box  in  place,  the  top  membrane  acted  as  a  gasket,  sealing  the  confining  fluid  around 
the  sides  from  the  fluid  on  the  top  of  the  specimen. 

5.2^  Instrumentation 

Two  types  of  instrumentation  were  employed  in  the  laboratory  investigation-  Pore 
water  pressure  measurements  were  made  using  transducers  attached  to  hypodermic  nee¬ 
dles.  This  allowed  for  direct  measurement  of  the  pore  water  pressure  within  the  interi¬ 
or  of  the  sample  while  causing  minimal  disturbance  to  the  specimen.  Small  ±  5  g 
accelerometers  were  placed  near  the  top  of  the  sample  and  at  the  sample  base  to  record 
table  (input)  motion  as  well  as  the  response  at  the  top  of  the  sample.  A  typical 
instrumentation  configuration  is  shown  in  Figure  37. 

5.2&  Test  Material 

Uniform  Ottawa  sand  was  used  in  all  tests  conducted  during  the  experimental 
program.  The  minimum  and  maximum  densities  were  determined  to  be  14.03  KN/cc 
and  15.99  KN/cc  respectively.  The  grain  size  distribution  for  the  test  sand  is  shown  in 
Figure  38. 
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Figure  38:  Ottawa  Sand  Gradation 
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5-3  TESTING  PROCEDURES 

5.3.1  Sample  saturation 

It  is  important  for  the  proper  measure  of  pore-pressure  rise  as  a  function  of 
shaking  that  the  sample  be  completely  saturated.  The  procedure  chosen  for  saturating 
the  test  specimens  was  to  boil  the  sample  in  deaired  water  under  a  vacuum  for  not 
less  than  30  minutes.  Each  sample  was  constructed  in  the  test  chamber  by  pluviation 
through  water.  This  method  was  chosen  because  it  has  been  our  experience  that  prepar¬ 
ing  specimens  this  way  will  consistently  yield  saturated,  uniform  samples  at  a  desired 
density.  In  order  to  increase  the  liquefaction  potential,  the  samples  tested  in  the  OSU 
studies  were  deposited  in  a  loose  state,  with  relative  densities  ranging  from  23%  to 
54%.  To  be  sure  full  saturation  was  achieved,  Skempton's  pore  water  pressure  parameter 
(  B  )  was  determined  prior  to  testing.  According  to  Black  and  Lee  [22],  values  of  B 
equal  to  1.0  signify  full  saturation  of  samples  of  low  relative  density.  A  value  for  B 
greater  than  0.95  signifies  an  acceptable  degree  of  saturation  for  sand  samples  having 
high  degrees  of  relative  density.  Considering  the  relative  densities  used  in  this  study, 
a  B  value  greater  than  or  equal  to  0.98  was  considered  indicative  of  full  saturation. 

5.3.2  Sample  Confinement 

Since  the  top  of  the  specimen  could  be  pressurized  independently  of  the  sides,  any 
ratio  of  vertical  to  horizontal  stress  could  have  been  used.  Two  different  different  rat¬ 
ios  were  tested.  The  response  of  samples  consolidated  under  isotropic  conditions  were 
the  majority  of  the  tests  performed  in  order  to  make  comparisons  with  other  experi¬ 
mental  data  practicaL  Several  samples  were  liquefied  after  being  consolidated  anisotropi- 
cally  (K  -  0.6). 


533  Input  Motions 

Two  types  of  input  motions  were  used  in  this  study.  The  first,  a  harmonic  input 
of  10  Hz  was  chosen  to  allow  for  comparison  between  published  liquefaction  data 
which  are  predominantly  periodic.  The  second,  a  random  amplitude  acceleration  time 
history  was  included  in  the  program  to  more  closely  simulate  the  type  of  motion  actu¬ 
ally  experienced  during  a  seismic  or  a  blast  induced  disturbance. 

53.4  Harmonic  motion 

A  nominal  10  Hz  table  motion  was  used  in  all  harmonic  tests.  The  10  Hz  fre¬ 
quency  was  chosen  because  it  was  high  enough  that  significant  acceleration  levels 
(>  2  g)  could  be  achieved  within  the  stroke  range  of  the  actuator  and  yet  was  well 
below  the  natural  frequency  of  the  unstrained  sample,  thus  permitting  an  assumption 
of  uniform  accelerations  throughout  the  height  of  the  sample.  As  can  be  seen  in  Figure 
39,  acceleration  time  histories  recorded  at  the  top  of  the  specimen  confirmed  that  the 
ratio  of  sample  to  table  accelerations  was,  in  fact,  approximately  equal  to  1.0  with  a 
phase  shift  consistent  with  a  wave  speed  of  about  1000  in/sec. 

De  Alba  [38]  had  constructed  their  sand  samples  using  a  reaction  mass  which  had 
been  placed  on  top  of  the  sample.  In  an  attempt  to  explain  experimentally  the  differ¬ 
ence  between  De  Alba's  and  our  results,  our  test  conditions  as  described  above  were 
modified.  The  modification  consisted  of  a  reaction  mass  similar  to  that  used  by  DeAlba 
being  added  to  the  test  specimen. 


5.3.5  Random  Motion 

A  standard  method  of  illustrating  the  frequency  content  of  an  acceleration  time 
history  is  by  means  of  a  Fourier  amplitude  spectrum.  The  Fourier  spectrum  is  an  indi¬ 
cation  of  the  final  energy  in  the  excitation  as  a  function  of  frequency.  In  the  context 
of  the  liquefaction  tests,  the  peaks  of  the  spectrum  represent  frequencies  at  which  rela¬ 
tively  large  amounts  of  energy  were  supplied  to  the  shaking  table/sand  system.  The 
Fourier  amplitude  spectrum  for  the  input  excitation  is  presented  in  Figure  40.  The 
spectra  show  the  dominant  frequencies  of  the  pink  noise  input  to  be  between  6  to  15 
Hz.  This  frequency  range  was  selected  to  symmetrically  bracket  the  harmonic  liquefac¬ 
tion  potential  data.  The  input  motions  were  derived  from  data  generated  on  an  IBM/PC 
XT  microcomputer  and  show  an  essentially  constant  amplitude  in  the  desired  frequency 
range.  Three  different  time  histories,  each  with  essentially  the  same  spectral  content, 
were  used  as  inputs  to  the  shaking  table. 

5-3-6  RESULTS 

The  studies  presented  in  this  chapter  were  conducted  to  provide  a  better  under¬ 
standing  of  the  liquefaction  phenomenon  through  large  scale  liquefaction  potential  test¬ 
ing,  and  to  provide  data  from  carefully  controlled  experiments  which  would  be  suitable 
for  use  in  model  verification.  Herein,  a  summary  covering  principal  findings  is  given. 
Details  are  available  in  items  1.11,  1.16,  1.18,  4.3  and  4.4  of  Appendix  B. 

The  results  of  the  harmonic  input  liquefaction  tests  performed  an  isotropically 
consolidated  samples  during  this  study  are  presented  in  Figure  41.  The  data  are  given 
in  terms  of  number  of  cycles  of  shaking  required  for  the  sample  reach  the  liquefied 
state  versus  the  cyclic  stress  level.  This  type  of  presentation  has  been  used  by  other 
investigators  [88,96,107,145],  and  is  presented  here  to  facilitate  a  comparison  between 
data  collected  in  this  study  with  those  presented  by  DeAlba  [38].  In  Figure  41,  it  can 
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Fourier  Spectrum  of  Random  Amplitude  Acceleration 
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be  seen  the  samples  tested  on  The  Ohio  State  University  facility  were  consistently  less 
susceptible  to  liquefaction  and  predictions  based  upon  our  results  would  result  in  much 
less  conservative  estimates  of  liquefaction  potential. 

When  the  results  of  the  anisotropic  tests  are  presented  in  the  same  format  as  the 
isotropic  test  data,  using  the  mean  normal  stress  to  represent  the  effective  confinement, 
the  curve  defining  the  number  of  cycles  to  liquefaction  plots  very  close  to  the  curve 
of  the  isotropic  results  as  can  be  seen  in  Figure  42. 

As  shown  in  Figure  43,  little  effect  on  the  liquefaction  potential  was  seen  exper¬ 
imentally  when  a  reaction  mass  similar  to  the  type  used  by  DeAlba  etaL  was  placed 
on  top  of  the  sample.  This  result  would  indicate  that  the  presence  or  absence  of  a 
reaction  mass  (regardless  of  its  height)  had  little  effect  on  the  liquefaction  potential  of 
the  sand  sample. 

The  results  of  the  different  tests  conducted  using  random  motions  are  presented 
in  Figure  44.  It  can  be  seen  from  this  figure  that  although  then  is  general  qualitative 
agreement  regarding  the  observed  liquefaction  potential  among  the  different  types  of 
random  dynamic  loadings  applied,  precise  quantitative  agreement  was  not  observed. 


Note:  Numbere  in  parentheses  indicate  items  of 
Appendix  B  which  describe  these  tests 


Liquefaction  Test  Results  for  Anisotropically  Consolidated  Samples 


Note:  Numbers  in  parentheses  indicate  items  of 
Appendix  B  which  describe  these  tests 


Figure  43:  Liquefaction  Test  Results  Obtained  with  a  Reaction  Mass  on  the  Sample 


Numbers  in  parentheses  indicate  items  of 
Appendix  B  which  describe  these  tests 


Liquefaction  Test  Results  for  Random  Amplitude  Excitation 


5.4  DISCUSSION 


The  motivation  for  the  test  program  conducted  at  The  Ohio  State  University  has 
been  the  generation  of  data  to  be  used  in  the  development  and  verification  of  mathe¬ 
matical  models.  This  purpose  was  different  from  any  of  the  previous  laboratory  pro¬ 
grams  in  that  the  chief  goal  of  those  programs  was  to  develop  empirical  relationships 
which  would,  for  a  limited  range  of  field  conditions,  allow  seismic  designs  to  be  made. 
Therefore,  in  the  current  study,  although  the  sample  boundary  conditions  applied  were 
designed  to  simulate  stresses  in  the  field  as  closely  as  was  practical,  it  was  imperative 
that  the  design  of  the  test  apparatus  allow  for  boundary  conditions  which  could  be 
clearly  defined,  that  the  material  constants  be  thoroughly  described,  and  the  density  of 
the  sample  be  as  uniform  as  practically  possible. 

Over  the  duration  of  the  project,  several  experimental  programs  have  been  con¬ 
ducted.  The  results  of  these  programs  have  shown  that  the  test  apparatus  used  was 
capable  of  operating  successfully  and  providing  consistent  results.  The  relatively  small 
offsets  between  the  points  identifying  the  onset  of  liquefaction  as  plotted  in  Figure  41 
can  be  explained  by  observing  that  similar  samples  were  often  of  slightly  different  rel¬ 
ative  densities.  It  was  also  shown  experimentally  that  tests  conducted  with  the  same 
testing  apparatus  continued  to  provide  consistent  results  even  with  the  addition  of  a 
reaction  mass  of  the  type  reported  by  DeAlba  [38}.  Therefore,  the  boundary  conditions 
imposed  on  the  top  of  the  sand  sample  by  this  reaction  mass  cannot  be  the  explanation 
as  to  why  the  testing  apparatus  used  in  this  study  gave  much  different  results  for 
liquefaction  potential  when  compared  with  the  results  presented  by  De  Alba  [38}  Other 
boundary  conditions  must  have  existed  in  their  samples  which  we  have  not  been  able 
to  reproduce  with  the  apparatus  used  in  this  study.  It  is  likely  that  the  precise  meth¬ 
od  used  by  DeAlba  [38]  to  attach  the  steel  shot  filled  bag  to  the  top  of  the  sample 
chamber  caused  localized  stress  concentrations  at  the  sand-shot  interface.  These  small 
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zones  of  high  stress  intensity  would  then  be  expected  to  generate  locally  high  pore 
pressures  resulting  in  premature  liquefaction. 

Although  some  soil  structure  rearrangement  is  certainly  likely  during  cyclic  load¬ 
ing,  we  find  it  difficult  to  attribute  the  eightfold  strength  increases  observed  by  Seed 
[158]  solely  to  the  change  in  structure  which  might  result  from  the  shaking  of  a 
homogeneous  sample,  particularly  while  maintaining  the  same  density.  What  is  more 
likely  is  that  the  pre-liquefaction  shaking  and  subsequent  drainage  of  a  specimen 
resulted  in  the  formation  of  a  more  uniform  sample.  Wolfe  [172]  has  shown  in  earlier 
tests  conducted  on  cubical  specimens  which  were  loaded  cyclically  to  liquefaction  and 
then  allowed  to  drain,  that  a  much  smaller  strength  increase  is  observed,  one  that  is 
consistent  with  the  observed  increase  in  relative  density  which  follows  shaking  and 
subsequent  reconsolidation. 

Three  different  band  limited  white  (or  pink)  noise  excitations  provided  the  shear 
stress  on  a  number  of  isotropically  consolidated  samples.  From  published  reports  avail¬ 
able  to  the  authors,  it  is  believed  that  these  were  the  first  such  tests  conducted  on 
large  scale  samples.  All  acceleration  time  histories  contained  the  same  dominant  frequen¬ 
cies  (6  to  15  Hz).  If,  however,  the  time  to  liquefaction  is  plotted  versus  maximum 
shear  maximum  shear  stress  ratio,  different  curves  for  the  different  input  motions 
emerge.  Therefore,  it  appears  that  the  use  of  the  maximum  shear  stress  within  the  soil 
sample  during  cyclic  loading  is  not  by  itself  an  appropriate  method  with  which  to 
characterize  the  effects  of  ground  shaking  on  liquefaction  potential. 

Several  potential  causes  of  the  differences  between  the  liquefaction  data  obtained 
in  this  study  and  the  data  reported  by  DeAlba  [38]  can  be  identified.  It  is  apparent 
that  liquefaction  is  essentially  a  local  phenomenon,  i-e.  local  irregularities  greatly  affect 
the  initiation  of  liquefaction.  Therefore,  non-uniformity  in  the  material  itself,  or  in 
the  load  applied  should  be  investigated. 
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5.4.1  Non-uniformity  of  sand  samples 

An  experimental  program  conducted  by  Mulilis  [107]  showed  the  importance  of 
the  method  of  sample  formation  on  the  liquefaction  resistance  of  a  specimen.  The  sam¬ 
ples  prepared  by  DeAlba  [38]  were  made  by  pluviation  through  air  and  subsequent  sat¬ 
uration,  whereas  all  the  specimens  prepared  in  the  present  study  were  saturated  by 
boiling  under  a  vacuum  then  deposited  entirely  under  water  at  a  constant  drop  height. 
The  data  presented  by  Mulilis  [107]  clearly  show  that  samples  made  by  these  two 
methods  demonstrate  different  cyclic  strengths.  They  attributed  the  strength  differences 
observed  to  significantly  different  soil  structures.  A  review  of  their  findings  shows 
that  although  the  wet  pluviation  technique  we  used  can  be  expected  to  result  in  high¬ 
er  cyclic  strengths  than  does  dry  pluviation,  it  is  not  likely  that  the  structural  differ¬ 
ences  obtained  alone  can  explain  the  magnitude  of  the  difference  we  have  observed 
between  the  two  sets  of  results.  Castro  [32]  attributes  the  cyclic  behavior  of  sands  to 
changes  in  the  void  spaces  at  the  local  level.  Marcuson  [95]  in  reporting  the  results  of 
an  experimental  program  designed  specifically  to  study  the  effects  of  sample  uniformity 
on  liquefaction  resistance,  observed  markedly  increased  cyclic  strengths  in  highly  uni¬ 
form  samples.  Furthermore,  Gilbert  [56]  has  recently  shown  that  samples  prepared  by 
pluviation  through  water  are  more  uniform  than  samples  prepared  by  the  popular 
method  of  moist  tamping. 

5.4.2  Nonunifonnities  due  to  testing 

Wood  [173]  observed  that  pore  pressures  measured  at  the  ends  of  the  triaxial  or 
simple  shear  sample  are  unlikely  to  be  useful  unless  tests  are  done  slowly  enough  to 
guarantee  pore  pressure  equalization  throughout  the  sample.  In  all  tests  conducted  dur¬ 
ing  this  study,  measurements  of  pore  pressure  rise  were  made  within  the  sample  itself 
not  at  a  sample  boundary.  Furthermore  pore  pressures  were  measured  at  more  than  one 


location  in  the  sample  during  every  test.  A  comparison  of  the  pore  pressure  time  histo¬ 
ries  during  the  test  indicates  that  at  any  instant  in  time  a  uniform  pressure 

throughout  the  central  portion  of  the  sample.  As  shown  in  Figure  44,  this  time 
dependent  but  spatially  independent  pore  pressure  was  typically  observed  throughout 
the  duration  of  the  test. 

5.4.3  Membrane  penetration 

Lade  [87]  and  Hernandez  studied  the  effects  of  membrane  penetration  on  the 
undrained  response  of  sands  in  triaxial  compression.  They  observed  that  the  existence  of 
membrane  penetration  in  samples  being  subjected  to  cyclic  loading  results  in  in  increase 

in  sample  volume  as  the  pore  pressures  in  the  sample  rise  during  shaking.  This  volume 

change,  if  not  accounted  for,  would  result  in  an  overestimation  of  the  liquefaction 
resistance  since  the  test  was  in  reality  a  partially  drained  test.  Lade  and  Hernandez 
cite  data  which  indicate  that  membrane  peneteration  is  negligible  for  sands  with  grain 
sizes  below  0.1  to  03.  mm,  or  for  effective  stresses  below  1.0  kg/sq  cm.  The  grain  size 
distribution  for  the  Ottawa  sand  used  in  this  study  was  given  in  Figure  38.  This  fig¬ 
ure  shows  that  the  sand  used  roughly  corresponds  to  sands  falling  in  the  range  where 
the  membrane  effects  would  be  condidered  negligible.  Also  effective  stresses  on  the 
sample  were  typically  on  the  order  of  1.0  kg/sqxm.  Nevertheless,  due  to  the  large  sam¬ 
ple  surface  covered  by  the  membrane,  it  was  felt  to  be  important  for  the  amount  of 
membrane  penetration  to  be  measured  accurately  in  order  to  assess  what  effect,  if  any, 
membrane  penetration  had  on  our  liquefaction  measurements. 

The  penetration  of  the  membrane  into  the  sample  was  measured  after  the  sample 
was  pressurized  to  the  test  conditions  and  before  the  internal  supporting  walls  were 
retracted.  The  test  chamber  was  attached  to  a  volume  change  measuring  device.  With 
the  isotropic  confining  pressure  constant  on  the  saturated  sand  sample,  the  back  pressure 
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was  incrementally  raised  to  simulate  generation  of  an  increase  in  pore  water  pressure 
during  actual  dynamic  testing.  Membrane  penetration  could  be  measured  as  movement 
of  the  membrane  out  from  between  the  sand  grains  and  recorded  in  the  form  of  vol¬ 
ume  change  on  a  manometer.  The  amount  of  membrane  penetration  measured  was 
consistently  less  than  1.75XKT4  percent  of  the  sample  volume.  This  amount  of  penetra¬ 
tion  was  considered  negligible. 

5.4.4  Nonuniformity  of  the  confining  pressure  at  the  sample  boundaries 

Since  in  the  shaking  table  configuration  used,  the  confining  water  must  be  accel¬ 
erated  along  with  the  sample,  there  was  some  concern  that  the  mass  of  the  water, 
which  was  arbitrarily  chosen,  could  be  affecting  the  distribution  of  pressures  on  the 
vertical  faces  of  the  sample.  Early  measurements  of  the  pressure  in  the  confining  fluid 
had  failed  to  show  any  variation  during  shaking,  but  it  was  felt  that  additional  test¬ 
ing  was  warranted.  Therefore,  in  a  specific  attempt  to  determine  experimentally  the 
effect  of  the  volume  of  confining  water  on  the  uniformity  of  the  confining  pressure 
and  therefore  on  liquefaction  potential,  plates  were  inserted  into  the  space  occupied  by 
the  confining  fluid.  The  effect  was  to  reduce  the  volume  of  water  by  more  than  50%. 
No  effect  on  liquefaction  potential  was  observed. 

The  results  of  a  test  program  in  which  saturated  sand  specimens  were  subjected 
to  harmonic  as  well  as  random  time  histories  of  base  accelerations  have  been  presented. 
The  data  show  that  loose  sands  can  be  liquefied  in  the  laboratory  and  that  the  resis¬ 
tance  to  liquefaction  is  strongly  dependent  upon  the  sample  boundary  conditions.  The 
test  apparatus  employed  in  the  program  described  herein  minimized  stress  irregularities 
at  sample  boundaries  and  should  therefore  be  seen  to  be  an  improvement  over  other 
laboratory  methods  for  determining  liquefaction  potential. 


Section  VI 


DISCUSSION 

6.1  OBJECTIVE  OF  THE  RESEARCH  PROGRAM 

The  objective  of  the  research  program  was  to  critically  examine  the  theoretical 
basis  of  the  equations  governing  the  behavior  of  saturated  soils  under  dynamic  loading 
in  order  to  identify/develop  appropriate  theory  which  would  properly  allow  for  soil- 
water  interaction  and  be  thermodynamically  consistent.  The  work  would  involve 
implementation  of  the  theory  or  theories  in  appropriate  solution  procedures.  A  program 
of  laboratory  investigation  was  included  to  constitute  a  reliable  database  for  verification 
of  the  theoretical  findings. 

6.2  ACCOMPLISHMENTS 

6.2.1  Review  of  Theories 

The  existing  theories  have  been  carefully  reviewed.  Limitations  of  the  theories 
and  their  differences  have  been  carefully  listed  and  discussed.  The  main  findings  of  the 
research  program  have  been  discussed  at  length  in  the  text  of  this  report  and  in  the 
publications  listed  in  Appendix  B.  A  brief  summary  of  the  conclusions  is  given  below. 

The  commonly  used  "engineering  approach"  to  study  of  liquefaction,  though  it  has 
been  successfully  utilized  for  the  study  of  many  case  histories,  cannot  be  directly 
extended  to  multi-dimensional  situations.  Even  for  simple  one-dimensional  cases,  it 
requires  considerable  "judgement"  on  the  part  of  the  engineer  in  addition  to  tedious 
Laboratory  investigations  on  the  dynamic  behavior  of  the  soil. 
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Biot's  theory  uses  some  ad  hoc  assumptions.  These  include  the  existence  of  an 
energy  function  for  the  mixture  and  the  existence  of  an  inertial  coupling  introduced 
through  a  kinetic  energy  expression  involving  product  of  the  velocities  of  the  soil  and 
the  pore-water.  There  are  questions  regarding  the  existence  and  the  nature  of  this 
coupling.  Definitely,  under  certain  circumstances  (e^»,  small  pore  size)  a  portion  of  the 
fluid  could  be  moving  effectively  with  the  soil.  This  coupling  would  only  constitute  a 
different  partitioning  of  the  total  mass  into  one  that  moves  with  the  soil  velocity  and 
the  remainder  which  moves  relaive  to  the  soil.  However,  the  introduction  of  a  "coupled 
mass”  appears  to  be  entirely  artificial  and  without  any  physical  basis.  Theories  of 
mixtures  have  been  developed  starting  from  different  assumptions.  Truesdell  assumed 
the  additivity  of  the  total  energy  of  the  constituents  and  proposed  artificial  definitions 
for  the  stresses  to  obtain  an  identity  of  form  between  the  equations  of  balance  for  the 
mixture  as  a  whole  and  those  for  the  individual  constituents.  It  is  apparent  that  the 
mixture  does  not  have  an  existence  as  a  continuum  in  motion  and,  therefore,  the  ques¬ 
tion  of  writing  equations  of  motion  for  it  should  not  arise  much  less  the  effort  to 
give  them  the  same  form  as  the  equations  for  each  constituent.  Truesdell's  "third  postu¬ 
late”  which  states  that  the  motion  of  the  mixture  is  governed  by  the  same  equations 
as  is  a  single  body,  is  unnecessary  and  irrelevant.  Green's  theory,  on  the  other  hand, 
assumes  the  additivity  of  stresses  and  fluxes.  However,  Green  as  well  as  Crochet,  intro¬ 
duced  energy  functions  for  the  mixture.  Their  work  does  not  necessarily  require  the 
notion  of  a  mixture  as  a  continuum  in  motion.  Bowen's  contention  that  Green's  theory 
is  a  special  case  of  Truesdell's  for  vanishing  relative  velocities  is  clearly  incorrect 
because  Green's  theory  is  no  less  general  than  Truesdell's.  There  are  differences  in  the 
meanings  attached  to  different  terms  even  though  the  form  of  the  equations  is  the 
same.  It  appears  that  balance  equations  ought  to  be  written  for  each  constituent  allow¬ 
ing  for  interaction.  Introduction  of  thermodynamic  quantities  associated  with  a  mix- 
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ture  as  a  continuum  in  motion  is  incorrect.  These  quantities,  in  relation  to  the  mixture, 
must  arise  as  a  consequence  of  the  quantities  associated  with  individual  constituents 
and  need  not  have  any  physical  interpretation. 

Truesdell  and  some  others  wrote  the  equations  of  balance  directly  considering  the 
motion  of  the  constituents  of  the  mixture  across  an  elementary  fixed  volume  in  space. 
The  contents  of  this  fixed  volume  change  constantly.  For  this  reason,  some  investigators 
have  proposed  writing  constitutive  equations  for  porosity.  This  would  be  unnecessary  if 
the  balance  laws  are  written  for  the  same  set  of  soil  particles  as  in  Gibson's  theory  of 
nonlinear  one-dimensional  consolidation. 

Constitutive  relations  are  required  for  the  diffusive  resistance  or  the  interaction 
force.  Truesdell ‘s  theory  of  mechanical  diffusion  includes  other  theories  as  specializa¬ 
tions.  Green's  theory  is  quite  similar. 

There  is  considerable  confusion  regarding  the  definition  of  partial  stresses.  Terza- 
ghi's  dual  definition  for  effective  stress  cannot  be  accepted.  The  effective  stress  is  dif¬ 
ferent  from  partial  stress  in  the  soil.  The  latter  includes  dependence  upon  the  kinemat¬ 
ics  of  the  fluid  in  addition  to  the  dependence  upon  the  kinematics  of  the  soil  whereas 

the  effective  stress  is  defined  as  the  part  of  the  soil  stress  which  depends  directly 

upon  the  strain  in  the  soil  skeleton.  Biot's  and  Green's  assumption  of  the  existence  of 

an  energy  function  for  the  mixture  would  lead,  in  the  linear  case,  to  a  symmetrical 

constitutive  coupling  between  the  soil  and  the  water.  However,  because  this  assumption 
is  questionable,  this  coupling,  even  if  it  exists,  need  not  be  symmetrical.  This  is  in  line 


with  Morland’s  views. 
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6J2J2  A  Dynamical  Theory  of  Saturated  Soil*. 

A  theory  of  dynamics  of  saturated  soils  was  developed.  This  theory  is  based  upon 
the  study  of  balance  of  a  fixed  set  of  particles  of  the  soil  contained  in  a  reference 
volume  in  the  initial  configuration.  This  reference  volume  moves  during  the  process 
and  changes  shape  as  well  The  description  of  deformation  and  motion  was  accom¬ 
plished  through  the  use  of  convected  coordinates  as  introduced  by  Novozhilov  [112]. 
This  eliminated  the  need  for  writing  constitutive  equations  for  porosity.  It  also  enabled 
a  clearer  definition  of  the  stresses  acting  on  the  representative  volume.  The  theory  may 
be  regarded  as  an  extension  of  Gibson's  theory  of  one-dimensional  consolidation  to  three- 
dimensions  and  inclusion  of  inertia  effects. 

6JL3  Development  of  Solution  Procedures. 

Truesdell's  as  well  as  Green's  theories,  for  the  case  of  small  motions  coincide  with 
Biot's  for  dynamics  of  saturated  soils.  The  new  dynamical  theory  of  saturated  soils 
would  also  coincide  with  Green's  theory  if  the  representative  volume  containing  the 
fixed  set  of  particles  undergoes  extremely  small  deformations.  Very  few  solutions  were 
available  for  the  problem.  Exact  solutions  had  been  developed  by  Biot  and  some  other 
investigators  for  some  simple  problems.  Numerical  solutions  had  been  attempted  but 
they  had,  in  general,  not  been  verified  against  exact  solutions  and  were  suspect.  In  the 
present  research  effort,  in  order  to  study  the  effectiveness  of  various  theories  in  mod¬ 
elling  dynamic  response  of  soil  systems,  it  was  necessary  to  systematically  develop 
solution  procedures.  These  included  exact,  semi-discrete,  as  well  as  numerical  solutions  to 
Biot's  formulation  of  the  problem  of  wave  propagation  in  saturated  soils  and  numerical 
solution  of  Seed's  and  Finn's  theory.  Solutions  were  also  developed  for  the  case  of  non¬ 


linear  material  behavior. 
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Exact  solutions  were  developed  for  several  cases  of  loading  of  an  infinite  column 
and  of  a  finite  layer  of  saturated  soil  For  the  infinite  column,  these  consisted  of 
integration  of  Garg's  fundamental  solution.  For  the  finite  layer  solutions  were  devel¬ 
oped  for  the  case  of  independent  specification  of  velocities  of  the  fluid  and  the  soil  at 
the  boundary.  In  case  of  a  sudden  application  of  velocity,  the  singularity  was  separated 
from  the  smooth  diffusion  and  the  two  solutions  superposed  for  the  linear  theory.  It 
was  noticed  that  Garg  had  made  some  assumptions  for  the  case  of  "weak  coupling” 
and  some  others  for  the  case  of  "strong  coupling”.  The  validity  of  these  assumptions 
has  been  carefully  examined  and  documented. 

Semi-discrete  solutions  were  developed  using  the  Laplace  Transform  technique  in 
conjunction  with  a  finite  element  discretization.  Effect  of  refinement  of  mesh  upon 
the  accuracy  of  the  results  was  examined.  The  eigenvalue  problem  becomes  extremely 
large  with  mesh  refinement.  This  approach,  though  useful  perhaps  for  benchmarking 
numerical  time-domain  solution  procedures,  was  seen  to  be  computationally  too  expen¬ 
sive. 


In  order  to  develop  finite  element  solution  procedures  in  a  systematic  manner,  the 
equations  of  motion  were  written  in  self-adjoint  form  in  an  appropriate  space.  Varia¬ 
tional  formulations  along  with  extensions  and  several  interesting  specializations  were 
developed.  Numerical  solutions  using  spatial  discretization  by  finite  elements  and  numer¬ 
ical  integration  over  the  time-domain  using  interpolation  schemes  of  suitable  order,  were 
developed  and  verified  against  some  exact  solutions.  Several  types  of  elements  were 
used  to  obtain  an  optimal  combination  of  accuracy  and  computational  economy.  For 
suddenly  applied  dynamic  disturbances,  it  was  found  necessary  to  use  "singularity"  ele¬ 
ments  to  properly  reproduce  the  wave  propagation. 
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For  Biot's  theory,  the  usual  form  is  the  two-field  formulation  in  which  the 
displacements  of  the  constituents  are  chosen  to  be  the  field  variables.  A  three-field  for¬ 
mulation  including  the  pore-water  pressures  as  the  additional  field  variable  was  written 
and  finite  element  solution  procedures  developed  for  the  same.  It  was  found  that  this 
scheme  gave  practically  the  same  results  as  the  two  field  formulation  for  the  problems 
in  which  the  velocities  were  specified.  However,  the  two-field  formulation  applies  the 
specified  fluid  pressure  condition  (e.g.  free-draining  boundary)  in  an  indirect  manner 
In  the  three  field  formulation  this  specification  is  direct. 

Computer  programs  for  two-dimensional  wave  programs  were  developed  but  could 
not  be  checked  against  exact  solutions  because  no  solutions  for  wave  propagation  in 
two-  or  three-  dimensional  saturated  soil  systems  were  available. 

Wave  propagation  through  nonlinear  saturated  soils  was  modelled  in  a  finite  ele¬ 
ment  based  computer  program.  No  exact  solutions  for  wave  propagation  in  nonlinear 
saturated  soils  are  available.  The  code  was  checked  against  two  solution  for  wave  pro¬ 
pagation  in  a  single  nonlinear  material  continuum. 

6.14  Laboratory  Investigations. 

Shaking  table  tests  were  conducted  on  saturated  samples  of  a  uniform  Ottawa 
sand.  The  use  of  a  shaking '  table  combined  with  the  large  size  of  the  samples  mini¬ 
mized  the  effects  of  stress  irregularities  at  sample  boundaries  and  the  method  used  to 
construct  the  samples  minimized  material  non-uniformities.  In  the  program  described  in 
this  report  both  harmonic  and  random  amplitude  table  acceleration  time  histories  were 
utilized.  A  significant  increase  in  the  resistance  to  liquefaction,  as  compared  with  the 
results  of  other  published  experimental  programs,  was  measured.  Potential  causes  for  the 
differences  in  test  results  were  investigated.  It  was  concluded  that  the  test  apparatus 
and  the  methods  employed  in  the  current  program  should  be  9een  as  an  improvement 
over  other  experimental  methods  for  determining  liquefaction  potential. 
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63  FUTURE  WORK 

Additional  work  on  research  in  the  general  area  of  dynamics  of  saturated  soils  is 
needed.  A  major  shortcoming  of  the  available  models  is  that  they  have  not  been  care¬ 
fully  tested  against  actual  experiments  in  the  field  or  in  the  laboratory.  The  reason  for 
this  is  that  experiments  simulating  simple  boundary  conditions  are  difficult  to  set-up 
and  most  experiments  that  are  convenient  to  carry  out  constitute  very  difficult 
boundary-value  problems  for  which  verified  computer  codes  are  difficult  to  come  by. 
There  is  need  to  design  test  set-ups  which  are  essentially  one-dimensional  wave  propa¬ 
gation  experiments  so  that  the  carefully  verified  solution  procedures  developed  in  the 
present  research  can  be  used  to  substantiate  the  theoretical  concepts  regarding  behavior 
of  saturated  soils.  On  the  other  hand,  there  is  need  to  develop  analytical  solutions  to 
two-  and  three-  dimensional  linear  as  well  as  nonlinear  problems  so  that  the  data 
from  shake-table  tests  which  are  essentially  two-  dimensional,  can  be  utilized  to 
verify  the  theoretical  models.  Any  of  the  models  would  need  material  properties  as 
input  data.  There  is  need  to  define  the  nonlinear  behavior  of  saturated  sands  very 
carefully  and  to  relate  it  to  the  propertries  of  the  single  materials  involved.  There 
has  been  difficulty  even  in  characterising  dry  sand  behavior  [llOj,  For  saturated 
materials,  the  role  of  the  pore-water  pressures,  the  existence  and  the  nature  of  the  con¬ 
stitutive  coupling  and  the  existence  or  otherwise  of  the  inertial  coupling  need  to  be 
studied.  The  relation  between  the  mechanics  of  the  particles  and  the  behavior  of  the 
soil  mass  needs  to  be  investigated.  The  soil  in  the  field  as  well  as  the  laboratory  is 
never  absolutely  uniform.  The  spatial  variation  in  material  properties  would  be  greater 
in  the  field  than  in  the  laboratory  where  the  conditions  can  be  controlled.  The  effect 
of  the  degree  of  randomness  on  the  overall  behavior  of  the  soil  is  extremely  important 
for  realistic  utilization  of  the  results  of  the  experiments  in  the  laboratory  for  the  situ¬ 
ation  in  the  field.  Also,  spatial  variation  in  material  properties  could  result  in  local 
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liquefaction  leading  to  perhaps  a  "domino”  effect  in  promoting  catastrophic  liquefaction 
in  the  soil  mass  in  the  field  even  though  the  soil  might  have  been  found  to  have 
high  resistance  to  liquefaction  in  laboratory  tests. 
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Appendix  A 

CONVECTED  COORDINATES 


Tin*  appendix  ««**«*«  a  aummary  of  aome  definition*,  relation*  and  formulae 
related  to  (he  uae  of  converted  coordinate*  in  the  mechanic*  of  continue.  The  discus¬ 
sion  follow*  Green  (65]  and  Fung  [47  ( 

A-l  NOTATION 

Indices.  which  may  either  be  Mbacnpta  or  superscripts.  such  as  x\  x,.  p1’.  etc. 

are  used  to  denote  component*  of  tenaora  of  various  order*.  A  single  element  $  hav¬ 

ing  no  indices  constitutes  *  system  of  tern  Older.  System*  of  element*  with  one  and 
two  indices,  ere  respectively,  termed  first  order  end  aeoood  order. 

The  summation  convention  used  throughout  this  report  implies  that  the  repetition 
of  an  index  (whether  superscript  or  subscript)  in  a  term  denotes  summation  with 

respect  to  that  index  over  the  range  1.2,  3. 

AJ  COORDINATE  TRANSFORMATION 

Let  9  denote  a  set  of  independent  variables,  whose  differentials  are  d?.  The 

mapping  of  9  into  another  set  of  variables  9  by  any  arbitrary  single- valued  function 
of  the  form 

9  *  9  C01.  91.  9J)  (A.l) 
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specific*  *  transformation  of  coordinates.  Here  we  have  used  the  same  notation  for  the 
function  and  its  value.  The  arbitrary  functions  are  aswmed  to  poaea  derivatives  upto 
the  order  required.  The  inverse  transformation  is  assumed  to  be  single  valued  and 
written  as: 

8*.  5J)  (A *2) 

The  differentials  d?  and  d7  are  related  by 
dJT -Cjd^ 

df-C^  (AJ) 

Evidently. 

CiCi-C,Ci**', 

where  8^  is  Krooecker's  delta.  For  reversibility  at  transformation  it  is  sufficient  that 

[M] 


C  «  Cl  *•  0  (A-5) 

While  the  transforms turns  of  the  differentials  in  (AJ)  is  linear,  the  transformations  of 
variables  in  (A.l)  and  (AJ)  are  not  necessarily  linear.  The  coordinate  transformation 
with  the  properties  dsernhert  above,  along  with  the  condition  (AAX  are  called  the 
admisnble  transformations  (47}  If  C  m  positive  everywhere,  then  a  right-handed  art  of 
coordinates  is  transformed  into  another  right-handed  set  and  the  transformation  is  prop¬ 
er.  In  this  work,  the  transformations  are  aswmed  to  be  admissible  and  proper,  unleei 


otherwise  stated. 
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A-3  CONTRA  VARIANT  AND  COVARIANT  VECTORS 

A  system  of  order  zero  has  only  a  single  component  £  in  the  variables  9*  and  a 
single  component  $  in  the  variables  9.  6  i»  scalar  if  $  ■  $  for  all  fr. 


Let  a  system  of  order  one  have  components  A'  and  X1.  in  the  variables 
9  and  9.  respectively.  Then,  if 

X,-C‘IA,-4?.A)  (Aj6) 

‘  & 

the  functions  X1  and  A1  art  contra vinant  components  of  a  tensor  of  order  one.  A  sys¬ 
tem  is  a  covariant  tensor  of  order  ooe  if  the  components  B,  m  the  variables  9  and  B, 
in  the  variables  9  are  related  by. 


B  -C°B  -  J^B 
'  1  jfiT  » 


(A.7) 


Contra  variant  tensor  components  are  indicated  by  superscripts  and  covmhant  components 
by  subscripts.  Transformations  of  contra  variant  and  covariant  components  of  second 
order  tensor*  are  given  by. 

(A.8) 


and 


B  *  C*  C°B 

ij  i  j  •• 


CA.9) 


(A-3)  shows  that  the  differentials  transform  according  to  the  law  of  contra  variant  ten- 
son.  The  use  of  upper  index  ts  appropriate  in  that  case.  The  variables  9  themselves 
are,  in  general,  neither  contra  variant  nor  covariant  and  the  position  of  their  index  is 
only  a  matter  of  convenience. 


A-4  CURVILINEAR  COORDINATES 


Consider  the  coordinates  x1  or  x,  referred  to  a  right  handed  orthogonal  system. 
Figure  45,  of  axes  in  a  three-dimensional  Euclidean  apace.  Let  an  admissible  transforma¬ 
tion  of  the  types  discussed  earlier  to  9  coordinates  be  defined  by 

d'«0'(xl.  xr  Xj)  (A.10) 

for  each  set  of  values  of  x*  there  exists  a  unique  set  of  values  of  9,  it  is  possi¬ 
ble  to  represent  the  Euclidean  space  by  the  variables  9,  instead  of  cartesian  system  x, 
The  relations 

9Cxr  xr  Xj) ■  constant  (All) 

with  i  •  1  to  3,  represent  three  families  of  coordinate  surfaces  and  the  point  of  inter¬ 
section  determines  a  point  P  in  space.  The  conditions  imposed  on  9  ensure  that  such  a 
point  is  uniquely  defined.  The  intersections  of  the  coordinate  surfaces  are  coordinate 
curves  and  9  are  curvilinear  coordinates. 

A-3  BASE  VECTORS  AND  METRIC  TENSORS 

Let  R  be  the  position  vector  of  a  point  P,  whose  coordinates  are  x,  or  x1  and  let 
dR  denote  the  infinitesimal  vector  PQ,  where  Q  has  the  coordinates  x,  +■  dx,  or  x‘  +  di' 
Then  (Figure  46) 

R*x.  e'  “x'e,  (A12) 

and 

dR  =  dx  e'  *  dx'  e.  (A13) 

i  i 

where  e,  and  e1  are  unit  vectors.  If  da  -  idRI  is  the  length  of  vector  PQ,  then 

(ds)2  =  dR .  dR  =  dx(  dx(  (A.14) 

If  0,,  0*  are  one  to  one  continuous  mappings  for  the  x,,  x1  systems. 


03  -curve 


-surface 
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d*1  «  &L  iji  ,  dx' «  4L  dd3 

&  & 

**>” 

Consequently 

dl  ■  dx  V  ■  &  e'd« -^d® 

$x'  1  J 

and  also 

dR  » di'  e  ■  -&-e.d0J*G1d&) 

1  &  '  1 

where 


(A1S) 

ex-lb) 


CA.17) 


(A.18) 


0»si«.  (A.I9) 

#  ‘ 

The  Kronecker  delta  and  the  permutation  eymbois  are  defined  in  rectangular 
cartesian  coordinate  system  as 


G  =G  .G 

>J  1  ) 


=  =  0  for  i*j  (  =1  for  i  =  j) 

(AJO) 

f 

0  (when  two  indices  are  equal) 

1  (when  i,  j,  k  are  in  even  permutation) 

|  —1  (when  i,  i  k  are  in  odd  permutation) 

(AJ1) 

in  general  coordinates  0*  are  as  follows: 

-&L&1 

(A-22) 

e*1  & 

_  afl* 

(A.23) 

3*  dx 

G‘  =G  G‘  =  ii-  =  s' 
J  j  ax“  ] 


(AJ4) 
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•i*  "  I-®-*  *  «,*  >/G 

V  'A 


|4®_I 


ax*  Jo 


(A-25) 


(A*26) 


when  O  it  the  determinant  tGJ  and  is  positive  for  any  proper  coordinate  system.  The 


vectora  O,  and  G1  are  connected  by  the  relations 
O'-0"O,  ;  O.-O,, O' 


(AJ7) 


Gu  is  the  Euclidean  metric  tensor  of  the  coordinate  system  and  Gu  is  the  associated 
metric  tensor  [47].  The  magnitudes  of  these  base  vectors  are. 


ig‘i  ■  Vg'.g1 


(AJ8) 


where  the  index  is  not  summed.  The  line  element  ds  in  (A.  14)  may  now  be  written 
in  the  form 


ds2  =  dR.dR  =  G  dO'd^^G^dfl  d0  =d0  d^ 

ij  i  j  • 


KJ>  PHYSICAL  MEANING  OF  COVARIANT  AND  CONTRA  VARIANT 
COMPONENTS  OF  A  VECTOR 

A  vector  w  in  terms  of  covariant  and  contra  variant  base  vectors  is 


w  =  w  G(  =  w(  G 


(A-29) 


(A30) 


where  by  (A.27), 


w.  =G  w*  w'  =  GiJ  w 
•  ij  j 


(AJ1) 
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According  to  (A-30),  if  w  is  represented  by  line  vectors,  then  w1  are  the  components 
of  w  in  the  direction  of  the  covariant  base  vectors  G,,  while  w,  are  the  component*  of 

w  in  the  direction  of  the  contravariant  base  vectors  O'  A  two-dimensional  illustration 
is  shown  in  Figure  47.  Some  other  results  are 

TW-GiJv‘wi»OlJviwJ»vw‘«v'  w. 

Iwl  »  v/w .  w  ■  V  Wj  w1  (A32) 

(AJO)  can  also  be  written  as 


w'VG^O,  w.^G' 

w  « - *  — ■ — — = —  (no  sum) 

w1  7G"B  are  components  of  w  resolved  in  the  direction  of  the  unit  vectors 


(A33) 

G, 

Tg; 


which  are  tangent  to  the  coordinate  curves.  A  similar  interpretation  for  w,  Vg“  is  pos¬ 
sible.  These  components  viz.  w1  -/G^ and  w,  Vg“  0  not  summed)  do  not  obey  the 
tensor  transformation  laws  and  are  not  components  of  tensors. 
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A.7  PARTIAL  DERIVATIVES  IN  CARTESIAN  COORDINATES 

Consider  1C,  a  comra variant  tenaor  of  order  one,  expressed  as  a  function  of  coor¬ 
dinate  system  The  relationahip  with  the  description  in  the  coordinate  system  x,  is 
given  by  CAj6)  ac 

X’Cx,.  lr  ij)«Ak(xlt  xr  (A34) 

d* 

Differentiating 


a?  ^axk  a?  *xk 

cajsj 

d*“d*‘  ax1  d?  a» 


If  the  two  cartesian  coordinate  systems  are  related  by  an  affine  transformation, 
i,  *  ay  x,  +  b^  au  and  b,  being  constants,  J£-  .  0  and  *  a^  .  Hence 

a*a» 


(AJ6) 

d?  a*‘  &  d*" 

(A-36)  implies  that  for  cartesian  systems,  related  by  an  affine  transformation,  the  par¬ 
tial  derivative  ]£!  follows  the  transformation  law  for  a  mixed  tenaor  of  order  two. 

d*‘ 

The  second  order  term  in  (A_35)  does  not  vanish  in  curvilinear  coordinates  and  thus 


the  partial  derivative  of  the  tenaor  doea  not  yield  a  higher  order  tensor. 
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AJ  COVARIANT  DIFFERENTIATION 

If  1  is  taken  to  be  a  scalar  function  of  the  general  coordinates  viz. 

R*R(0,.  d3)  (AJ7) 

then 


The  differentiation  of  base  vectors  G,  yields 
G  =  3*  =  R  =Ra 

*0  Qgi  A)  *jl 

Recalling  (A  19)  and  noting  that  e,  form  a  set  of  constant  base  vectors, 
G  ■ 

J  ae  V  r 

which,  upon  use  of  (A.19),  becomes 


(A38) 


(AJ9) 


(A.40) 


where 


and 


g  =  r  G*  =  rr  g 

y  *j»  ‘J  r 


.2  r  __ r 

r  =  v  *-  &- 

'*  a^a^  a®* 


rV0"1* 


(A-22)  and  (A.42)  give 


(A.41) 


(A.42) 


(A.43) 


(A.44) 


Similarly,  (A-29)  implies 

G|j  =  -rijrGr  (A.45) 

The  symbols  riJ(  and  rTti  are  Christoff  el  symbols  of  the  first  and  second  kind,  respec¬ 
tively.  If  w  is  a  vector 
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ffcL  =  &!L  =  Cr  =  Cr 

a?  d0r  e®1  1 5®r 


1  .r 


(A.46) 


so  that  w^  transform  according  to  the  covariant  transformation.  Also,  from  CA-30), 


w.  =  w^Gr  +  wrGr. 


=  wr,Gr  +  wrGrj 


Using  (A.41)  and  (A.4S) 


w.  =  wrl.  G  =  w  I.Gr 

A  i  r  r  i 


where 


wi.  =  w  .  —  r* .  w 

r  i  rj  rl  ■ 


(A47) 

(A-48) 

(A.49) 

(A^O) 


The  expressions  w'l,  and  wrl,  are  the  covariant  derivatives  of  wr  and  respectively,  of 
the  vector  w.  Since  wi  transform  according  to  covariant  transformation,  it  follows 
from  (A.48)  that  the  covariant  derivatives  of  a  vector  form  a  tensor  of  order  two 
Similarly,  we  can  write  expressions  for  derivatives  of  second  order  tensors  in  the  fol¬ 
lowing  form. 


ijr 


+  r‘ 


rm 


A 

A 


-r“  a. 

mj  jr  lm 

“j  +  rj  Ain 

rm 


(A^l) 


A.9  GEOMETRIC  INTERPRETATION  OF  COVARIANT  DERIVATIVES 

Consider  a  vector  field  w  associated  with  every  point  in  space  in  the  region 
under  consideration.  Let  the  vector  w  at  a  point  P  (0\  O2,  8s)  be 

w(P)  =  wi(P)Gi(P)  (A-52) 

At  a  neighboring  point  Q  (d1  +dd\  Q2  +d92,  d3  -Kid5),  the  vector  is 
w  (Q)  =  w  (P)  +d  w  (P) 
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=  [w‘(P)  +  dw‘(P)]  [G(  (P)  +  dG,  (P)  ] 

Taking  the  limit  as  d0*  -*  0,  we  get, 

dw  =  [w‘  +  dw1]  [Gj  +  dG,]  —  w1  Gj  =  w1  dG.  +  dw*  G, 
and  the  derivative 


(A.S3) 


(A.S4) 


(A-55) 


Thus,  the  derivative  of  w  consists  of  two  parts  viz.  one  due  to  the  variation  of  the 
components  of  w*  as  the  coordinates  0“  are  changed  and  the  other  arising  from  the 
change  of  the  base  vectors  G,  as  the  position  of  the  point  0*  is  changed.  Use  of  (A.41) 
in  (A-55)  gives 

=  +  w“  r'.G 

J 


=[iZL  +  w‘  rr.i  g 

&  r 


=  wrI.G 


(A.56) 


Thus  the  contravariant  derivative  w'lj  represents  the  components  of  referred  to  the 
base  vectors  Gr. 

A.10  GRADIENT,  DIVERGENCE  AND  CURL  IN  CURVILINEAR  COORDINATES 
The  gradient  (grad  <f>),  divergence  (div  F)  and  curl  (curl  A)  where 
^  is  a  scalar  and  F,  A  are  vectors  are  given  by 

grad  ^  =  <f>  r  Gr  (A.S7) 

div  F  =  Frl  (A.58) 


curl  A  =  €  Atlf  Gt 


(A.59) 
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A.U  KINEMATICS  AND  KINETICS 
A.11.1  Introductory 

The  particles  of  matter  that  occupy  regions  of  Euclidean  space  form  a  body.  A 
given  body,  B,  may  occupy  different  regions  at  different  times.  For  convenience,  we 
denote  the  configuration  at  time  r  =  0,  as  a  reference  configuration  Q.  The  fired 
region  of  space  C,  is  occupied  by  particles  of  the  body.  Now,  if  this  material  moves 
so  that  at  a  subsequent  time  r-t,  it  occupies  a  new  region  of  space  C,  then  the  body 
is  said  to  have  undergone  a  motion.  Mathematically,  this  can  be  represented  by  a  series 
of  coordinate  transformations.  In  particular,  the  proper  transformations  assumed  in  this 
work  ensure  that  the  axioms  of  continuity  and  impenetrability  are  followed.  As  the 
continuum  moves  from  one  configuration  to  another,  the  matter  in  the  neighborhood  of 
each  point  is  translated  and  rotated  as  a  rigid  body  and  is  strained.  Strain  of  an  ele¬ 
mental  volume  is  that  part  of  the  relative  motion  between  neighboring  particles  that  is 
not  due  to  the  rigid  motion.  In  this  section  various  measures  of  strain  due  to  relative 
displacements  are  described  with  reference  to  a  single  material. 

A1U  Geometric  Relations 

Consider  a  continuous  three-dimensional  body  in  the  reference  configuration  C0  at 
time  r  =  0.  Let  a  system  of  coordinates  a,  be  so  chosen  that  a  point  P  of  the  body  is 
described  by  P  (a,,  a*  a*).  At  time  r  =  t,  let  the  body  be  in  configuration  C,  having 
undergone  a  motion.  The  particle  P,  originally  in  C,  now  be  at  Q  in  C,  with  coordi¬ 
nates  (x,,  xT  Xj).  The  coordinate  systems  a,  and  x,  may  be  curvilinear  and  need  not  be 
coincident.  They  both  describe  a  Euclidean  space.  Figure  48. 

The  admissible  transformations  of  the  type  described  in  Section  A.2  are 
x.  =  x.(al,  a2,  a3)  (A.67) 


which  has  a  unique  inverse 


18 


ai=ai(ll-  XX*  X3} 


(A.68) 


Consider  an  infinitaaimal  line  element  PP  in  the  reference  configuration  C,  such 


that  P  is  given  by  P (a,  +  da,). Then,  ds,,  the  length  of  PP,  is 


ds2  =  a11da'daJ 


(A  .69) 


Here  a^  is  the  Euclidean  metric  for  the  coordinate  system  a,.  Let  a  corresponding  line 


segment  in  the  configuration  C  be  QQ  and  its  length  ds.  Then 
ds2  »Gjjdxidxj 


(A.70) 


where  G(J  is  the  Euclidean  metric  tensor  for  the  coordinate  system  z,.  Recalling  (A.15), 


The  change  in  length  is: 

ds2  —  ds2  =  [G  Jdaidaj 

°  “  9a1  $aJ  ^ 

or  equivalently, 

ds2  —  ds2  =[G.  —  a  ^.J^l]dxidri 

o  1  >j  »n  i  j 

ax  a* 

Defining  strain  tensors 

E..  =  I[G  ^^-a.J 

lJ  2  3a‘  $aJ  ^ 

e.  =  — [G  .  —  a  2L-2L1 
,J  2  ,J  mn  ax1  & 


(A.71) 


CA.72) 


(A.73) 


(A.74) 


we  get 


ds2  —  ds2  =  2  ESj  da'  daj 


ds2  —  ds2  =  2  e,  dx1  dr* 

O  i| 


(A.75) 


(A.76) 


E,.  and  e„,  are,  respectively.  Green's  strain  tensoT  and  Almansi  strain  tensor  [47]. 
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A.11.3  Description  of  Deformation 

Generally,  the  initial  and  the  deformed  configurations  are  described  with  reference 
to  a  fixed  rectangular  cartesian  frame. 

However,  an  alternative  introduced  by  Novozhilov  [112]  is  to  use  a  frame  of  ref¬ 
erence  which  moves  continuously  with  the  body  and  deforms  in  such  a  way  that  the 
numerical  values  of  the  position  coordinates  of  the  particles  of  the  body  remain  the 
same  throughout  the  motion.  This  is  convenient  for  finite  deformation  [47]. 

For  fluid  motion,  the  interest  is  centered  not  on  the  displacements  of  the  particles 
or  their  velocities  but  the  velocity  distribution  in  the  volume  occupied  by  the  fluid. 
The  most  convenient  mathematical  framework,  for  this  purpose  is  the  Eulerian  coordi¬ 
nates  [112]. 


Consider  a  rectangular  cartesian  coordinate  system  a,  shown  in  Figure  49.  The 
position  vector  of  a  typical  point  P#  in  the  reference  configuration  Q  is 

r  =  a.  e. .  (A-77) 

where  e,  are  unit  vectors  along  fixed  axes.  Let  a  typical  point  P0  in  C0  occupy  posi¬ 
tion  P  in  the  new  configuration  C  at  time  r  =  (O^r^t).  The  position  vector  of  P 
referred  to  the  same  fixed  axes  is 

E  =  z.ei  (A.78) 

Components  of  the  displacement  vector  are: 

ui  =  r  -  a.  (A.79) 

Also,  the  one  to  one  correspondence  between  points  in  C  and  Q  implies 

Z-SZ-U,.  a2,  a3)  (A. 80) 


(A.81) 
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where  a,  and  z,  are  admissible  transformations.  Here  we  do  not  differentiate  between  a 
function  and  its  value.  The  difference  will  be  clear  in  the  context. 

The  velocity  ▼  of  a  particle  referred  to  the  fixed  system  of  axes,  is 

y  =  i  =  Vie.  (A.82) 

It  follows  from  (A.79)  that 

v.  =  z.  =  u.  (A.83) 

Similarly,  the  acceleration  f  of  a  particle  referred  to  the  fixed  system  of  axes  is 

f  =  f.e.  =  v  =  i  (A.84) 

and 

f.  =  z.  ss'uj  =  Vj  (A.85) 

The  initial  state  in  C„  may  also  be  described  by  a  general  curvilinear  set  of  coor¬ 
dinates  x,  so  that, 

a^ajCi,,  x2,  x3,  r  =  0)  (A.86) 

where  a,  is  single-valued  and  continuously  differentiable  as  many  times  as  required. 
We  may  imagine  this  curvilinear  system  x,  to  move  continuously  with  the  body  as  we 
pass  from  the  reference  state  C,  at  r  =0  to  the  state  C  at  r  =  t.  The  triplet  of  real 
numbers  x,  are  merely  the  labels  that  we  assign  to  the  positions  corresponding  to  the 
material  particles  of  the  body  in  the  reference  configuration  C„  The  values  of  x, 
remain  unchanged  for  the  new  position  P  in  C.  The  coordinates  x,  may  be  rectangular 
or  curvilinear  in  but  are  curvilinear,  in  general,  in  C.  Thus  the  coordinates  of  P 
with  reference  to  the  invariant  convected  coordinates  are 

Zj  =  z.(x,,  x2,  x3,  r  =  t)  (A.87) 

Using  (A.15)  and  (A.86),  we  may  define  a  contra  variant  vector  dx1,  in  Q,  as 

dx'  =  daj  da'  =  dxj  (A.88) 


From  (A.80)  and  (A.88) 


Noting  that  the  convected  frame  x,  is  rectangular  in  C„,  the  base  vectors  g,,  g1  for  the 
system  x,  in  the  body  at  C#  are 


Ki  =  r.i  =  c, 

(A.94) 

gi.gj  =  e'.ej  =  S‘ 

(A.95) 

The  metric  tensors  gtJ,  g1J  are 

8tj  =  *i  •  8j  =  giJ  =  g‘.gj  =  S,j 

(A.96) 

Similarly,  the  base  vectors  G, ,  G1  and  metric  tensors  Gu ,  GIJ  for  the  curvilinear  system 
x,  in  the  body  at  time  t  are 


G  =  R  G.Gj  =  8j 

I  ,1  1  i 

G  =  G  .G  =  &-  3*L 

,J  '  J  ^  ^ 

G‘j  =  G' .  GJ  =  -3L 

8zr  & 


(A.97) 

(A.98) 

(A.99) 


The  strain  tensor  in  convected  coordinate  system  x,  is  defined  as 
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2  «  2 

For  line  elements  corresponding  to  the  vectors  r  and  R  viz.  ds,  and  ds, respectively. 


(A.100) 


dso  =  «iJdxidlJ  =  8ijdxidxJ 

ds2  =  G,jdxidxJ 


CA.101) 

(A.102) 


and 


ds2  —  ds2  =  2  ydx’  dr* 

O  9  lj 


(A.103) 

yu  may  be  expressed  in  terms  of  the  displacement  vector  u  or  its  components 
with  respect  to  the  base  vectors  g,  (Le.  e,)  or  Gr  We  have 

Gi  =  R,ia  r.i +  =  ei +  “i  =  +  umP  e»  =  (8»i +  U»P  e»  (A-104) 

Hence, 

VU  =  T ' u. j +  *j •  u.i +  u.i •  ®.  P ^  I (Gi • tt. j +  Gi •  u.i ~ u.i •  “P  (A*105) 


The  displacement  vector  u  may  be  expressed  in  terms  of  base  vectors  of  C,  and  C. 
Thus, 


u  =  ue 

m  m 


u  =  u  e 

.1  ,  i  m 


and 


u  =  U  Gm  =  UmG 

m  to 


u  =  U  I  Gm  =  Uml.  G 

,1  mi  i  m 


(A.106)  and  (A.105)  give 

y  =  4-  [u.  +u.  +  u  u  J 

Similarly,  (A.  107)  along  with  (A.97)  through  (A.99)  gives 

y  =  i[U.I  +U  I.-U  I.  U  lj 
'ij  2  1  J  J  *  m  I  m  r 


(A.106) 


(A.107) 


(A.108) 


(A.109) 
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A.  11.4  Velocity  and  Acceleration  in  Conroe  ted  Coordinates 

When  referred  to  the  state  C,  the  expressions  for  velocity  and  acceleration  vec¬ 
tors  with  respect  to  x,  coordinates  are 

v  *  u*  Oj  ®  u(  e'  *  Uj  e(  CA-110) 

f  =  ulei  =  u[e,  =  ule(  (A-lll) 

The  vectors  ▼  and  f  can  also  be  resolved  in  terms  of  the  base  vectors  C, ,  G1  of  the 
configuration  C  Several  alternative  representations  are  possible.  Thus, 

▼  =  v"G  *®v  Gb  (A.112) 

is  m 


(A.113) 

(A.114) 

(A.11S) 

(A.116) 

(A-117) 

(A118) 

(4.119) 


(A.  120) 


2 


f  -  lii  -  u‘  (A.121) 

y 

A.Ui  Qu|«  is  Volume  end  Area 

The  convened  coordinate  system  is  rectangular  cartesian  in  the  initial  state  and 
curvilinear  in  the  deformed  state.  Consider  an  elementary  rectangular  perailelopiped 
enclosed  by  the  six  surfaces  viz. 

z  ■  constant  and  x,  +  dx)  ■  constant  (i-1  to  3)  (A.122) 

The  vectors  s,  dx,  (no  sum)  form  the  sides  of  the  parallepiped,  the  volume  of  which  is 
given  by. 

dV0-dxidlJdxj  -dx‘dxJdxJ  (A.123) 

After  motion  and  deformation,  the  rectangular  parallepiped  is  deformed  into  a  perailo- 
piped  with  sides  defined  by  G,  dx‘  (no  sum)  1  *  1, 2, 3_  The  new  volume  is 

dV-C|.(01xG,>dx1dxJdxJ  (A.124) 

but 

O1.(0JxGJ)i»lrJ»IGtJl»  VG  (A-125) 

Hence 

dV  »  yGdi,  dijdjj  ■  i/GdV0  (A.126) 

If  dAM  denotes  an  elementary  area  in  the  reference  configuration  on  the  coordinate 
plane  x,  *  constant,  then 

dAJ0 » '•2xe,l dZjdZj  » dx2 dZj  * dx2 dx3  (A.127) 

After  deformation,  the  originally  plane  area  dA„  becomes  a  curved  surface  dA,  with 
vectors  Gj  dxJ  and  G,  dx1  as  its  sides.  The  area 


dA,  *  G,xG,ldx2dx3 


(A.128) 


(A.129) 


W**®*-  ^e23tO‘=VGG‘ 

The  magnitude  IG2xG,l  is 

IG2xG3l  =  IVGG’I  =  VgG11  '  (A.130) 

Hence 

dA,  =  VGG11  dx2  dx3  CA.131) 

In  general, 

dA,*  VGGU dxJdxk  (no  sum  on  i  and  i^  ja*k)  (A.132) 

in  which  dA,  denote  the  areas  in  the  deformed  body  which  in  the  undeformed  state 
has  values  dA*  (i  -  1  to  3). 

A.  11.6  Kinetics 

This  discussion  follows  [65].  Let  the  surfaces  x,=  constant  at  a  point  P  in  the 

deformed  state  C  form  a  tetrahedron.  The  edges  of  the  tetrahedron  are  formed  by  the 

coordinate  curves  PP,  of  length  ds,  and  the  curves  P,Pj,  P^P,  P,?,,  (Figure  50).  Let  the 

surfaces  x,  *  constant  of  the  tetrahedron  have  the  areas  dA,.  These  may  be  represented 
vectorially  by 

,  dA.  , 

G*  — (A.133) 

Jo* 

Ql 

where  — — -  is  the  unit  base  vector.  Also,  the  area  of  P,PJ\  denoted  by  dA,  is  repre- 

Vg5 

rented  vectorially  by  ndA,  where  n  is  the  unit  normal  to  the  surface.  Hence,  since 
the  area  P,PjP,  is  vectorially  equivalent  to  the  areas  dA, 
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so  that  if  n,  are  the  covariant  components  with  respect  to  the  base  vectors  G',  it  fol¬ 
lows  that, 

i  G‘dA 

-0dA=^ 

or  n.  Vg**  dA  =  dA.  (A.  135) 


Let  t  denote  the  stress  per  unit  area  of  the  surface  at  a  point  P  in  the  deformed 


body.  Cauchy's  equation  of  motion,  following  Green  [65l  is 


■f tiA+{ 


p(F  — f)dV  =  0 


(A.136) 


where  V  an  the  arbitrary  volume  in  the  body  in  the  state  C  and  is  bounded  by  a 
closed  surface  A.  Also,  p  is  the  density  of  the  body  in  C  F  and  f  are  the  body  force 
and  acceleration  vectors,  respectively.  Applying  the  equation  of  motion  to  an  infinitesi¬ 
mal  tetrahedron  PP1PJPJ  we  have,  in  the  limit,  keeping  the  direction  of  n,  the  normal 


to  the  area,  fixed. 


t  dA  « t.  dA; 


(A.137) 


Note  that  dA,  are  the  areas  of  the  surfaces  of  the  tetrahedron  under  consideration  and 
t,  are  the  stress  vectors  associated  with  these  elemental  areas.  Volume  forces  and  iner¬ 
tia  terms  acting  on  the  tetrahedron  do  not  appear  in  this  equation  since  they  are  of 
higher  order  of  smallness  than  the  surface  forces.  (A.135)  and  (A.136)  give 


t  =  ZMi^ 


(A.138) 


The  stress  vector  t  is  invariant  under  transformation  of  coordinates  and  n,  are  compo¬ 
nents  of  covariant  vector.  It  implies  that  t,  Vg“  is  a  contra  variant  type  transformation 


t.  VGS  =  rijG.  =  TiGj 
•  j  j 


(A.139) 


where  r*  and  are  componenets  of  the  associated  contravariant  and  mixed  tensors  of 
second  order,  called  stress  vectors.  (A.138)  and  (A.139)  give 


t=f!£=T‘VrTSn>Gi 


(A.  140) 


where 


Tj  «  t.  Vgg“  =  \fG  riJ G.  =  n/Gt* Gj 


(A.141) 


T,  is  introduced  for  convenience.  An  element  of  area  at  a  point  x,  in  the  body  C  is 
V GGU dxJ dxk  (i  not  summed  and  i^j^k)  and  the  force  across  this  element  is, 

t.  VGGSdxjdxk  =  T,dxjdxk  (A.142) 


The  conditions  at  the  boundary  surface  of  the  body,  at  which  surface  tractions  are  pre¬ 
scribed,  require 


t  =  P  =  P'Gi  =  PjGJ 


(A.  140)  gives 

rV-pi 


(A.  143) 


(A.  144) 


t*  n  =  P 

i  i  i 


(A.145) 


It  is  generally  assumed  that  the  geometry  of  the  reference  configuration  is  known  a 
priori,  it  is  often  convenient  to  define  the  state  of  stress  at  a  point  by  relating  to  its 
position  in  the  initial  configuration.  The  stress  vector,  t,  is  referred  to  a  surface  S  at 
time  t  in  body  C  and  measured  per  unit  area  of  S.  The  stress  tensor  r*1  referred  to  x, 
coordinates  in  C  and  is  measured  per  unit  area  of  these  coordinate  surfaces. 


If  st,  is  a  stress  vector  across  the  x,  surface  in  C,  but  measured  per  unit  area  of 
the  corresponding  x,  surface  in  Cp  then  use  of  (A.142),  yields 
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If  is  the  stress  vector  measured  per  unit  area  of  the  undeformed  body,  associated 
with  a  surface  in  the  deformed  body,  whose  unit  normal  in  its  undeformed  position  is 

0t  =  sViGj  =  »>iGj  CA.147) 

where  „n  =  ^  g1  =  ,n‘  g,.  (A.  139)  gives 

0tiV7  =  0tjV?  =  sijGj  (A.148) 

(A.148)  and  (A.146)  give 

Tj  =  Qt.  V  g  gu  =  >/g  sij G.  =  slj  Gj  as  >/g  =  1  for  a  cartesian  frame  (A.149) 

where 

sij=>/GTij  (A.150) 

s*J  are  componets  of  the  second  Piola-Kirchoff  stress  tensor  [114].  (A.104)  along  with 

(A.141)  and  (A.149)  gives 

Tj  **  >/G rij (8mj  +  ujem  =  slj (8mj  +  eo  -  >/G P‘m ^ * S‘ffi CA.151) 
The  stresses  P*.  =  r11(fi,J  +  um>P  are  measured  per  unit  area  of  x,  surfaces  in  C  but  are 
referred  to  the  base  vectors  e,  in  Q.  The  stresses  S^s^CS^  +  u^p  are  measured  per 
unit  area  of  the  x,  surfaces  in  C,  and  are  referred  to  the  base  vectors  e,  in  C,  These 
are  not  symmetric.  The  stress  tensor  defined  by  components  is  called  the  first 
Piola-Kirchoff  stress  tensor. 
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